OSDN Git Service

These files are integrated into libtetrabz.c
authorMitsuaki Kawamura <kawamitsuaki@gmail.com>
Thu, 24 Mar 2022 11:54:49 +0000 (20:54 +0900)
committerMitsuaki Kawamura <kawamitsuaki@gmail.com>
Thu, 24 Mar 2022 11:54:49 +0000 (20:54 +0900)
12 files changed:
src/c/__init__.c [deleted file]
src/c/libtetrabz_common.c [deleted file]
src/c/libtetrabz_common.h [deleted file]
src/c/libtetrabz_dbldelta.c [deleted file]
src/c/libtetrabz_dblstep.c [deleted file]
src/c/libtetrabz_dos.c [deleted file]
src/c/libtetrabz_fermigr.c [deleted file]
src/c/libtetrabz_occ.c [deleted file]
src/c/libtetrabz_polcmplx.c [deleted file]
src/c/libtetrabz_polstat.c [deleted file]
src/c/test.c [deleted file]
src/c/test.py [deleted file]

diff --git a/src/c/__init__.c b/src/c/__init__.c
deleted file mode 100644 (file)
index 35566e0..0000000
+++ /dev/null
@@ -1,31 +0,0 @@
-//\r
-// Copyright (c) 2014 Mitsuaki Kawamura\r
-//\r
-// Permission is hereby granted, free of charge, to any person obtaining a\r
-// copy of this software and associated documentation files (the\r
-// "Software"), to deal in the Software without restriction, including\r
-// without limitation the rights to use, copy, modify, merge, publish,\r
-// distribute, sublicense, and/or sell copies of the Software, and to\r
-// permit persons to whom the Software is furnished to do so, subject to\r
-// the following conditions:\r
-//\r
-// The above copyright notice and this permission notice shall be included\r
-// in all copies or substantial portions of the Software.\r
-//\r
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS\r
-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.\r
-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY\r
-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,\r
-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE\r
-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
-//\r
-from .libtetrabz_dos import dos\r
-from .libtetrabz_dos import intdos\r
-from .libtetrabz_occ import occ\r
-from .libtetrabz_occ import fermieng\r
-from .libtetrabz_polstat import polstat\r
-from .libtetrabz_fermigr import fermigr\r
-from .libtetrabz_dbldelta import dbldelta\r
-from .libtetrabz_dblstep import dblstep\r
-from .libtetrabz_polcmplx import polcmplx\r
diff --git a/src/c/libtetrabz_common.c b/src/c/libtetrabz_common.c
deleted file mode 100644 (file)
index ec85a22..0000000
+++ /dev/null
@@ -1,594 +0,0 @@
-/*
-// Copyright (c) 2014 Mitsuaki Kawamura
-//
-// Permission is hereby granted, free of charge, to any person obtaining a
-// copy of this software and associated documentation files (the
-// "Software"), to deal in the Software without restriction, including
-// without limitation the rights to use, copy, modify, merge, publish,
-// distribute, sublicense, and/or sell copies of the Software, and to
-// permit persons to whom the Software is furnished to do so, subject to
-// the following conditions:
-// 
-// The above copyright notice and this permission notice shall be included
-// in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include <stdio.h>
-#include <stdlib.h>
-#include "libtetrabz.h"
-
-void libtetrabz_initialize(
-    int ng[3],
-    double bvec[3][3],
-    double wlsm[20][4],
-    int **ikv
-    )
-    {
-  /*
-  define shortest diagonal line & define type of tetragonal
-  */
-  int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
-  double bvec2[3][3], bvec3[4][3], bnorm[4];
-  double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
-                         {0.0, 1440.0, 0.0, 30.0},
-                         {30.0, 0.0, 1440.0, 0.0},
-                         {0.0, 30.0, 0.0, 1440.0} },
-  wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
-                 {-28.0, -38.0, 7.0, 17.0},
-                 {17.0, -28.0, -38.0, 7.0},
-                 {7.0, 17.0, -28.0, -38.0}},
-  wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
-                 {9.0, -56.0, 9.0, -46.0},
-                 {-46.0, 9.0, -56.0, 9.0},
-                 {9.0, -46.0, 9.0, -56.0}},
-  wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
-                 {7.0, -38.0, -28.0, 17.0},
-                 {17.0, 7.0, -38.0, -28.0},
-                 {-28.0, 17.0, 7.0, -38.0}},
-  wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
-                 {-18.0, -18.0, -18.0, 12.0},
-                 {12.0, -18.0, -18.0, -18.0},
-                 {-18.0, 12.0, -18.0, -18.0}};
-
-  for (i1 = 0; i1 < 3; i1++)
-    for (i2 = 0; i2 < 3; i2++)
-      bvec2[i1][i2] = bvec[i1][i2] / (double) ng[i1];
-
-  for (i1 = 0; i1 < 3; i1++) {
-    bvec3[0][i1] = -bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
-    bvec3[1][i1] = bvec2[0][i1] - bvec2[1][i1] + bvec2[2][i1];
-    bvec3[2][i1] = bvec2[0][i1] + bvec2[1][i1] - bvec2[2][i1];
-    bvec3[3][i1] = bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
-  }
-  /*
-  length of delta bvec
-  */
-  for (i1 = 0; i1 < 4; i1++) {
-    bnorm[i1] = 0.0;
-    for (i2 = 0; i2 < 3; i2++)
-      bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
-  }
-  itype = 0;
-  for (i1 = 1; i1 < 4; i1++)
-    if (bnorm[i1] < bnorm[itype]) itype = i1;
-  /*
-  start & last
-  */
-  for (i0 = 0; i0 < 4; i0++) {
-    ivvec0[i0] = 0;
-    for (i1 = 0; i1 < 4; i1++)divvec[i0][i1] = 0;
-    divvec[i0][i0] = 1;
-  }
-  ivvec0[itype] = 1;
-  divvec[itype][itype] = -1;
-  /*
-  Corners of tetrahedron
-  */
-  it = 0;
-  for (i0 = 0; i0 < 3; i0++) {
-    for (i1 = 0; i1 < 3; i1++) {
-      if (i1 == i0) continue;
-      for (i2 = 0; i2 < 3; i2++) {
-        if (i2 == i1 || i2 == i0) continue;
-
-        for (i3 = 0; i3 < 3; i3++) {
-          ivvec[it][0][i3] = ivvec0[i3];
-          ivvec[it][1][i3] = ivvec[it][0][i3] + divvec[i0][i3];
-          ivvec[it][2][i3] = ivvec[it][1][i3] + divvec[i1][i3];
-          ivvec[it][3][i3] = ivvec[it][2][i3] + divvec[i2][i3];
-        }
-
-        it += 1;
-      }
-    }
-  }
-  /*
-  Additional points
-  */
-  for (i1 = 0; i1 < 6; i1++) {
-    for (i2 = 0; i2 < 3; i2++) {
-      ivvec[i1][4][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][1][i2];
-      ivvec[i1][5][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][2][i2];
-      ivvec[i1][6][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][3][i2];
-      ivvec[i1][7][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][0][i2];
-
-      ivvec[i1][8][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][2][i2];
-      ivvec[i1][9][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][3][i2];
-      ivvec[i1][10][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][0][i2];
-      ivvec[i1][11][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][1][i2];
-
-      ivvec[i1][12][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][3][i2];
-      ivvec[i1][13][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][0][i2];
-      ivvec[i1][14][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][1][i2];
-      ivvec[i1][15][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][2][i2];
-
-      ivvec[i1][16][i2] = ivvec[i1][3][i2] - ivvec[i1][0][i2] + ivvec[i1][1][i2];
-      ivvec[i1][17][i2] = ivvec[i1][0][i2] - ivvec[i1][1][i2] + ivvec[i1][2][i2];
-      ivvec[i1][18][i2] = ivvec[i1][1][i2] - ivvec[i1][2][i2] + ivvec[i1][3][i2];
-      ivvec[i1][19][i2] = ivvec[i1][2][i2] - ivvec[i1][3][i2] + ivvec[i1][0][i2];
-    }
-  }
-
-  for (i1 = 0; i1 < 4; i1++) {
-    for (i2 = 0; i2 < 4; i2++) {
-      wlsm[i2][i1] = wlsm1[i1][i2] /= 1260.0;
-      wlsm[i2 + 4][i1] = wlsm2[i1][i2] /= 1260.0;
-      wlsm[i2 + 8][i1] = wlsm3[i1][i2] /= 1260.0;
-      wlsm[i2 + 12][i1] = wlsm4[i1][i2] /= 1260.0;
-      wlsm[i2 + 16][i1] = wlsm5[i1][i2] /= 1260.0;
-    }
-  }
-  /*
-  k-index for energy
-  */
-  nt = 0;
-  for (i2 = 0; i2 < ng[2]; i2++) {
-    for (i1 = 0; i1 < ng[1]; i1++) {
-      for (i0 = 0; i0 < ng[0]; i0++) {
-
-        for (it = 0; it < 6; it++) {
-
-          for (ii = 0; ii < 20; ii++) {
-            ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
-            ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
-            ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
-            for (i3 = 0; i3 < 3; i3++) if (ikv0[i3] < 0) ikv0[i3] += ng[i3];
-            ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
-          }
-          nt += 1;
-        }
-      }
-    }
-  }
-}
-
-
-void libtetrabz_tsmall_a1(
-    double *e,
-    double e0,
-    double *v,
-    double tsmall[4][4]
-) {
-  /*
-  Cut small tetrahedron A1
-  */
-  double a10, a20, a30;
-  a10 = (e0 - e[0]) / (e[1] - e[0]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-  a30 = (e0 - e[0]) / (e[3] - e[0]);
-
-  *v = a10 * a20 * a30;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 1.0 - a10;
-  tsmall[0][2] = 1.0 - a20;
-  tsmall[0][3] = 1.0 - a30;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = a10;
-  tsmall[1][2] = 0.0;
-  tsmall[1][3] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = a20;
-  tsmall[2][3] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = 0.0;
-  tsmall[3][3] = a30;
-}
-
-
-void libtetrabz_tsmall_b1(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron B1
-  :param e: numpy
-  :return:
-  */
-  double a13, a20, a30;
-  a13 = (e0 - e[3]) / (e[1] - e[3]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-  a30 = (e0 - e[0]) / (e[3] - e[0]);
-
-  *v = a20 * a30 * a13;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 1.0 - a20;
-  tsmall[0][2] = 1.0 - a30;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 0.0;
-  tsmall[1][2] = 0.0;
-  tsmall[1][3] = a13;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = a20;
-  tsmall[2][2] = 0.0;
-  tsmall[2][3] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = a30;
-  tsmall[3][3] = 1.0 - a13;
-}
-
-
-void libtetrabz_tsmall_b2(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron B2
-  :param e:
-  :return:
-  */
-  double a21, a31;
-  a21 = (e0 - e[1]) / (e[2] - e[1]);
-  a31 = (e0 - e[1]) / (e[3] - e[1]);
-
-  *v = a21 * a31;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 0.0;
-  tsmall[0][2] = 0.0;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 1.0;
-  tsmall[1][2] = 1.0 - a21;
-  tsmall[1][3] = 1.0 - a31;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = a21;
-  tsmall[2][3] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = 0.0;
-  tsmall[3][3] = a31;
-}
-
-void libtetrabz_tsmall_b3(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron B3
-  :param e:
-  :return:
-  */
-  double a12, a20, a31;
-  a12 = (e0 - e[2]) / (e[1] - e[2]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-  a31 = (e0 - e[1]) / (e[3] - e[1]);
-
-  *v = a12 * a20 * a31;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 1.0 - a20;
-  tsmall[0][2] = 0.0;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 0.0;
-  tsmall[1][2] = a12;
-  tsmall[1][3] = 1.0 - a31;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = a20;
-  tsmall[2][2] = 1.0 - a12;
-  tsmall[2][3] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = 0.0;
-  tsmall[3][3] = a31;
-}
-
-
-void libtetrabz_tsmall_c1(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron C1
-  :param e:
-  :return:
-  */
-  double a32;
-  a32 = (e0 - e[2]) / (e[3] - e[2]);
-
-  *v = a32;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 0.0;
-  tsmall[0][2] = 0.0;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 1.0;
-  tsmall[1][2] = 0.0;
-  tsmall[1][3] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = 1.0;
-  tsmall[2][3] = 1.0 - a32;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = 0.0;
-  tsmall[3][3] = a32;
-}
-
-
-void libtetrabz_tsmall_c2(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron C2
-  :param e:
-  :return:
-  */
-  double a23, a31;
-  a23 = (e0 - e[3]) / (e[2] - e[3]);
-  a31 = (e0 - e[1]) / (e[3] - e[1]);
-
-  *v = a23 * a31;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 0.0;
-  tsmall[0][2] = 0.0;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 1.0;
-  tsmall[1][2] = 1.0 - a31;
-  tsmall[1][3] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = 0.0;
-  tsmall[2][3] = a23;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = a31;
-  tsmall[3][3] = 1.0 - a23;
-}
-
-void libtetrabz_tsmall_c3(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][4]
-)
-{
-  /*
-  Cut small tetrahedron C3
-  :param e:
-  :return:
-  */
-  double a23, a13, a30;
-  a23 = (e0 - e[3]) / (e[2] - e[3]);
-  a13 = (e0 - e[3]) / (e[1] - e[3]);
-  a30 = (e0 - e[0]) / (e[3] - e[0]);
-
-  *v = a23 * a13 * a30;
-
-  tsmall[0][0] = 1.0;
-  tsmall[0][1] = 1.0 - a30;
-  tsmall[0][2] = 0.0;
-  tsmall[0][3] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 0.0;
-  tsmall[1][2] = a13;
-  tsmall[1][3] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = 0.0;
-  tsmall[2][3] = a23;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = a30;
-  tsmall[3][2] = 1.0 - a13;
-  tsmall[3][3] = 1.0 - a23;
-}
-
-
-void libtetrabz_triangle_a1(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][3]
-)
-{
-  /*
-  Cut triangle A1
-  :param e:
-  :return:
-  */
-  double a10, a20, a30;
-  a10 = (e0 - e[0]) / (e[1] - e[0]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-  a30 = (e0 - e[0]) / (e[3] - e[0]);
-
-  *v = 3.0 * a10 * a20 / (e[3] - e[0]);
-
-  tsmall[0][0] = 1.0 - a10;
-  tsmall[0][1] = 1.0 - a20;
-  tsmall[0][2] = 1.0 - a30;
-  tsmall[1][0] = a10;
-  tsmall[1][1] = 0.0;
-  tsmall[1][2] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = a20;
-  tsmall[2][2] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2] = a30;
-}
-
-void libtetrabz_triangle_b1(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][3]
-) {
-  /*
-  Cut triangle B1
-  :param e:
-  :return:
-  */
-  double a30, a13, a20;
-  a30 = (e0 - e[0]) / (e[3] - e[0]);
-  a13 = (e0 - e[3]) / (e[1] - e[3]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-
-  *v = 3.0 * a30 * a13 / (e[2] - e[0]);
-
-  tsmall[0][0] = 1.0 - a20;
-  tsmall[0][1] = 1.0 - a30;
-  tsmall[0][2] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = 0.0;
-  tsmall[1][2] = a13;
-  tsmall[2][0] = a20;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = a30;
-  tsmall[3][2] = 1.0 - a13;
-}
-
-
-void libtetrabz_triangle_b2(
-    double *e,
-    double e0,
-    double *v,
-    double tsmall[4][3]
-)
-{
-  /*
-
-  :param e:
-  :returns:
-      - x - first
-      - y - second
-  */
-  double a12, a31, a20;
-  a12 = (e0 - e[2]) / (e[1] - e[2]);
-  a31 = (e0 - e[1]) / (e[3] - e[1]);
-  a20 = (e0 - e[0]) / (e[2] - e[0]);
-
-  *v = 3.0 * a12 * a31 / (e[2] - e[0]);
-
-  tsmall[0][0] = 1.0 - a20;
-  tsmall[0][1] = 0.0;
-  tsmall[0][2] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = a12;
-  tsmall[1][2] = 1.0 - a31;
-  tsmall[2][0] = a20;
-  tsmall[2][1] = 1.0 - a12;
-  tsmall[2][2] = 0.0;
-  tsmall[3][0] = 0.0;
-  tsmall[3][1] = 0.0;
-  tsmall[3][2]= a31;
-}
-
-void libtetrabz_triangle_c1(
-  double *e,
-  double e0,
-  double *v,
-  double tsmall[4][3]
-)
-{
-/*
-Cut triangle C1
-*/
-  double a03, a13, a23;
-  a03 = (e0 - e[3]) / (e[0] - e[3]);
-  a13 = (e0 - e[3]) / (e[1] - e[3]);
-  a23 = (e0 - e[3]) / (e[2] - e[3]);
-
-  *v = 3.0 * a03 * a13 / (e[3] - e[2]);
-
-  tsmall[0][0] = a03;
-  tsmall[0][1] = 0.0;
-  tsmall[0][2] = 0.0;
-  tsmall[1][0] = 0.0;
-  tsmall[1][1] = a13;
-  tsmall[1][2] = 0.0;
-  tsmall[2][0] = 0.0;
-  tsmall[2][1] = 0.0;
-  tsmall[2][2] = a23;
-  tsmall[3][0] = 1.0 - a03;
-  tsmall[3][1] = 1.0 - a13;
-  tsmall[3][2] = 1.0 - a23;
-}
-
-void eig_sort(
-    int n, //!< [in] the number of components
-    double *key, //!< [in] Variables to be sorted [n].
-    int *swap //!< [out] Order of index (sorted)
-)
-{
-  int i, j, k, min_loc;
-  double min_val;
-
-  for (i = 0; i < n; ++i) swap[i] = i;
-
-  for (i = 0; i < n - 1; ++i) {
-    min_val = key[i];
-    min_loc = i;
-    for (j = i + 1; j < n; ++j) {
-      if (min_val > key[j]) {
-        min_val = key[j];
-        min_loc = j;
-      }
-    }
-    if (key[i] > min_val) {
-      /*
-       Swap
-      */
-      key[min_loc] = key[i];
-      key[i] = min_val;
-
-      k = swap[min_loc];
-      swap[min_loc] = swap[i];
-      swap[i] = k;
-    }
-  }/*for (i = 0; i < n - 1; ++i)*/
-}/*eig_sort*/
diff --git a/src/c/libtetrabz_common.h b/src/c/libtetrabz_common.h
deleted file mode 100644 (file)
index 0505522..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-/*
- Copyright (c) 2014 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-void libtetrabz_initialize(int ng[3], double bvec[3][3], double wlsm[20][4], int** ikv);
-void libtetrabz_tsmall_a1(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_b1(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_b2(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_b3(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_c1(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_c2(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_tsmall_c3(double *e, double e0, double *v, double tsmall[4][4]);
-void libtetrabz_triangle_a1(double *e, double e0, double *v, double tsmall[4][3]);
-void libtetrabz_triangle_b1(double *e, double e0, double *v, double tsmall[4][3]);
-void libtetrabz_triangle_b2(double *e, double e0, double *v, double tsmall[4][3]);
-void libtetrabz_triangle_c1(double *e, double e0, double *v, double tsmall[4][3]);
-void eig_sort(int n, double *key, int *swap);
diff --git a/src/c/libtetrabz_dbldelta.c b/src/c/libtetrabz_dbldelta.c
deleted file mode 100644 (file)
index 46cfa49..0000000
+++ /dev/null
@@ -1,266 +0,0 @@
-//
-// Copyright (c) 2014 Mitsuaki Kawamura
-//
-// Permission is hereby granted, free of charge, to any person obtaining a
-// copy of this software and associated documentation files (the
-// "Software"), to deal in the Software without restriction, including
-// without limitation the rights to use, copy, modify, merge, publish,
-// distribute, sublicense, and/or sell copies of the Software, and to
-// permit persons to whom the Software is furnished to do so, subject to
-// the following conditions:
-// 
-// The above copyright notice and this permission notice shall be included
-// in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-//
-#include "libtetrabz_common.h"
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-static void libtetrabz_dbldelta2(
-    int nb,
-    double **ej,
-    double **w
-) {
-  /*
-  2nd step of tetrahedron method.
-  */
-  int i3, ib, indx[3];
-  double a10, a20, a02, a12, v, e[3], e_abs;
-
-  for (ib = 0; ib < nb; ib++) {
-
-    e_abs = 0.0;
-    for (i3 = 0; i3 < 3; i3++) {
-      e[i3] = ej[ib][i3];
-      if (e_abs < fabs(e[i3])) e_abs = fabs(e[i3]);
-    }
-    eig_sort(3, e, indx);
-
-    if (e_abs < 1.0e-10) {
-      printf("Nesting ##\n");
-    }
-
-    if ((e[0] < 0.0 && 0.0 <= e[1]) || (e[0] <= 0.0 && 0.0 < e[1])) {
-
-      a10 = (0.0 - e[0]) / (e[1] - e[0]);
-      a20 = (0.0 - e[0]) / (e[2] - e[0]);
-
-      v = a10 / (e[2] - e[0]);
-
-      w[ib][indx[0]] = v * (2.0 - a10 - a20);
-      w[ib][indx[1]] = v * a10;
-      w[ib][indx[2]] = v * a20;
-    }
-    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
-
-      a02 = (0.0 - e[2]) / (e[0] - e[2]);
-      a12 = (0.0 - e[2]) / (e[1] - e[2]);
-
-      v = a12 / (e[2] - e[0]);
-
-      w[ib][indx[0]] = v * a02;
-      w[ib][indx[1]] = v * a12;
-      w[ib][indx[2]] = v * (2.0 - a02 - a12);
-    }
-    else {
-      for (i3 = 0; i3 < 3; i3++)
-        w[ib][i3] = 0.0;
-    }
-  }
-}
-
-void dbldelta(
-  int ng[3],
-  int nk,
-  int nb,
-  double bvec[3][3],
-  double** eig1,
-  double** eig2,
-  double*** wght
-)
-{
-  /*
-  Main SUBROUTINE for Delta(E1) * Delta(E2)
-  */
-  int it, ik, ib, i20, i4, j3, jb, ** ikv, indx[4];
-  double wlsm[20][4], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[4][3], thr;
-
-  thr = 1.0e-10;
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (ik = 0; ik < 6 * nk; ik++) {
-    ikv[ik] = ikv[0] + ik * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ej1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  ej1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    ei1[i4] = ei1[0] + i4 * nb;
-    ej1[i4] = ej1[0] + i4 * nb;
-  }
-
-  w1 = (double***)malloc(nb * sizeof(double**));
-  w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
-  w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
-    }
-  }
-
-  ej2 = (double**)malloc(nb * sizeof(double*));
-  ej2[0] = (double*)malloc(nb * 3 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    ej2[ib] = ej2[0] + ib * 3;
-  }
-
-  w2 = (double**)malloc(nb * sizeof(double*));
-  w2[0] = (double*)malloc(nb * 3 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w2[ib] = w2[0] + ib * 3;
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] = 0.0;
-
-  for (it = 0; it < 6 * nk; it++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ib = 0; ib < nb; ib++) {
-        ei1[i4][ib] = 0.0;
-        ej1[i4][ib] = 0.0;
-      }
-    for (i20 = 0; i20 < 20; i20++) {
-      for (i4 = 0; i4 < 4; i4++) {
-        for (ib = 0; ib < nb; ib++) {
-          ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
-          ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
-        }
-      }
-    }
-
-    for (ib = 0; ib < nb; ib++) {
-
-      for (i4 = 0; i4 < 4; i4++)
-        for (jb = 0; jb < nb; jb++)
-          w1[ib][i4][jb] = 0.0;
-
-      for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] < 0.0 && 0.0 <= e[1]) {
-
-        libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (jb = 0; jb < nb; jb++)
-            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
-          libtetrabz_dbldelta2(nb, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
-        }
-      }
-      else if (e[1] < 0.0 && 0.0 <= e[2]) {
-
-        libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (jb = 0; jb < nb; jb++)
-            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
-          libtetrabz_dbldelta2(nb, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
-        }
-
-        libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (jb = 0; jb < nb; jb++)
-            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
-          libtetrabz_dbldelta2(nb, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
-        }
-      }
-      else if (e[2] < 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (jb = 0; jb < nb; jb++)
-            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
-          libtetrabz_dbldelta2(nb, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j3 = 0; j3 < 3; j3++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
-        }
-      }
-      else {
-        continue;
-      }
-    }
-    for (i20 = 0; i20 < 20; i20++)
-      for (ib = 0; ib < nb; ib++)
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] /= (6.0 * (double)nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(ej1[0]);
-  free(ej1);
-  free(ej2[0]);
-  free(ej2);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1); 
-  free(w2[0]);
-  free(w2);
-}
diff --git a/src/c/libtetrabz_dblstep.c b/src/c/libtetrabz_dblstep.c
deleted file mode 100644 (file)
index 86be968..0000000
+++ /dev/null
@@ -1,376 +0,0 @@
-/*
- Copyright (c) 2014 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include "libtetrabz_common.h"
-#include <math.h>
-#include <stdlib.h>
-
-static void libtetrabz_dblstep2(
-    int nb,
-    double ei[4],
-    double **ej,
-    double **w1
-) {
-  /*
-  Tetrahedron method for theta( - de)
-  */
-  int i4, j4, ib, indx[4];
-  double v, thr, e[4], tsmall[4][4];
-
-  thr = 1.0e-8;
-
-  for (ib = 0; ib < nb; ib++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      w1[ib][i4] = 0.0;
-
-    for (i4 = 0; i4 < 4; i4++) e[i4] = - ei[i4] + ej[ib][i4];
-    eig_sort(4, e, indx);
-
-    if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
-      /*
-      Theta(0) = 0.5
-      */
-      v = 0.5;
-      for (i4 = 0; i4 < 4; i4++)
-        w1[ib][i4] += v * 0.25;
-    }
-    else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
-      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-    }
-    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
-
-      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-
-      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-
-      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-    }
-    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
-
-      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-
-      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-
-      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j4 = 0; j4 < 4; j4++)
-          w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
-    }
-    else if (e[3] <= 0.0) {
-      for (i4 = 0; i4 < 4; i4++)
-        w1[ib][i4] += 0.25;
-    }
-  }
-}
-
-void dblstep(
-    int ng[3],
-    int nk,
-    int nb,
-    double bvec[3][3],
-    double **eig1,
-    double **eig2,
-    double ***wght
-)
-{
-    /*
-    Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
-    */
-  int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4];
-  double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr;
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (ik = 0; ik < 6 * nk; ik++) {
-    ikv[ik] = ikv[0] + ik * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ej1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  ej1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    ei1[i4] = ei1[0] + i4 * nb;
-    ej1[i4] = ej1[0] + i4 * nb;
-  }
-
-  ej2 = (double**)malloc(nb * sizeof(double*));
-  ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    ej2[ib] = ej2[0] + ib * 4;
-  }
-
-  w1 = (double***)malloc(nb * sizeof(double**));
-  w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
-  w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
-    }
-  }
-
-  w2 = (double**)malloc(nb * sizeof(double*));
-  w2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w2[ib] = w2[0] + ib * 4;
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] = 0.0;
-
-  thr = 1.0e-8;
-
-  for(it = 0; it < 6*nk; it++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ib = 0; ib < nb; ib++) {
-        ei1[i4][ib] = 0.0;
-        ej1[i4][ib] = 0.0;
-      }
-    for (i20 = 0; i20 < 20; i20++) {
-      for (i4 = 0; i4 < 4; i4++) {
-        for (ib = 0; ib < nb; ib++) {
-          ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
-          ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
-        }
-      }
-    }
-    
-    for (ib = 0; ib < nb; ib++) {
-
-      for (i4 = 0; i4 < 4; i4++)
-        for (jb = 0; jb < nb; jb++)
-          w1[ib][i4][jb] = 0.0;
-
-      for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] <= 0.0 && 0.0 < e[1]) {
-
-        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[1] <= 0.0 && 0.0 < e[2]) {
-
-        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[2] <= 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_dblstep2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[3] <= 0.0) {
-        for (i4 = 0; i4 < 4; i4++) {
-          ei2[i4] = ei1[i4][ib];
-          for (jb = 0; jb < nb; jb++)
-            ej2[jb][i4] = ej1[i4][jb];
-        }
-        libtetrabz_dblstep2(nb, ei2, ej2, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            w1[ib][i4][jb] += w2[jb][i4];
-      }
-      else {
-        continue;
-      }
-    }
-    for (i20 = 0; i20 < 20; i20++)
-      for (ib = 0; ib < nb; ib++)
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] /= (6.0 * (double) nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(ej1[0]);
-  free(ej1);
-  free(ej2[0]);
-  free(ej2);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-  free(w2[0]);
-  free(w2);
-}
diff --git a/src/c/libtetrabz_dos.c b/src/c/libtetrabz_dos.c
deleted file mode 100644 (file)
index 6275998..0000000
+++ /dev/null
@@ -1,281 +0,0 @@
-/*
- Copyright (c) 2022 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include "libtetrabz_common.h"
-#include <stdlib.h>
-
-void dos(
-  int ng[3],
-  int nk,
-  int nb,
-  int ne,
-  double bvec[3][3],
-  double** eig,
-  double* e0,
-  double*** wght) {
-  /*
-  Compute Dos : Delta(E - E1)
-  */
-  int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
-  double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][3];
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (ik = 0; ik < 6 * nk; ik++) {
-    ikv[ik] = ikv[0] + ik * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (ii = 0; ii < 4; ii++) {
-    ei1[ii] = ei1[0] + ii * nb;
-  }
-
-  w1 = (double***)malloc(nb * sizeof(double**));
-  w1[0] = (double**)malloc(nb * ne * sizeof(double*));
-  w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * ne;
-    for (ie = 0; ie < ne; ie++) {
-      w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
-    }
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (ie = 0; ie < ne; ie++)
-        wght[ik][ib][ie] = 0.0;
-
-  for (it = 0; it < 6 * nk; it++) {
-
-    for (ii = 0; ii < 4; ii++)
-      for (ib = 0; ib < nb; ib++)
-        ei1[ii][ib] = 0.0;
-    for (jj = 0; jj < 20; jj++)
-      for (ii = 0; ii < 4; ii++)
-        for (ib = 0; ib < nb; ib++)
-          ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
-
-    for (ib = 0; ib < nb; ib++) {
-
-      for (ie = 0; ie < ne; ie++)
-        for (ii = 0; ii < 4; ii++)
-          w1[ib][ie][ii] = 0.0;
-
-      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
-      eig_sort(4, e, indx);
-
-      for (ie = 0; ie < ne; ie++) {
-
-        if ((e[0] < e0[ie] && e0[ie] <= e[1]) || (e[0] <= e0[ie] && e0[ie] < e[1])) {
-
-          libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 3; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
-
-        }
-        else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
-
-          libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 3; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
-
-          libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 3; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
-        }
-        else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
-
-          libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 3; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
-        }
-        else {
-          continue;
-        }
-      }
-    }
-    for (ii = 0; ii < 20; ii++)
-      for (ib = 0; ib < nb; ib++)
-        for (ie = 0; ie < ne; ie++)
-          for (jj = 0; jj < 4; jj++)
-            wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (ie = 0; ie < ne; ie++)
-        wght[ik][ib][ie] /= (6.0 * (double)nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-}
-
-void intdos(
-  int ng[3],
-  int nk,
-  int nb,
-  int ne,
-  double bvec[3][3],
-  double** eig,
-  double* e0,
-  double*** wght
-) {
-  /*
-  Compute integrated Dos : theta(E - E1)
-  */
-
-  int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
-  double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][4];
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (ik = 0; ik < 6 * nk; ik++) {
-    ikv[ik] = ikv[0] + ik * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (ii = 0; ii < 4; ii++) {
-    ei1[ii] = ei1[0] + ii * nb;
-  }
-
-  w1 = (double***)malloc(nb * sizeof(double**));
-  w1[0] = (double**)malloc(nb * ne * sizeof(double*));
-  w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * ne;
-    for (ie = 0; ie < ne; ie++) {
-      w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
-    }
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (ie = 0; ie < ne; ie++)
-        wght[ik][ib][ie] = 0.0;
-
-  for (it = 0; it < 6 * nk; it++) {
-
-    for (ii = 0; ii < 4; ii++)
-      for (ib = 0; ib < nb; ib++)
-        ei1[ii][ib] = 0.0;
-    for (jj = 0; jj < 20; jj++)
-      for (ii = 0; ii < 4; ii++)
-        for (ib = 0; ib < nb; ib++)
-          ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
-
-    for (ib = 0; ib < nb; ib++)
-      for (ie = 0; ie < ne; ie++)
-        for (ii = 0; ii < 4; ii++)
-          w1[ib][ie][ii] = 0.0;
-
-    for (ib = 0; ib < nb; ib++) {
-
-      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
-      eig_sort(4, e, indx);
-
-      for (ie = 0; ie < ne; ie++) {
-
-        if ((e[0] <= e0[ie] && e0[ie] < e[1]) || (e[0] < e0[ie] && e0[ie] <= e[1])) {
-          libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-        }
-        else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
-
-          libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-          libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-          libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        }
-        else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
-
-          libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-          libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-          libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
-          for (ii = 0; ii < 4; ii++)
-            for (jj = 0; jj < 4; jj++)
-              w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        }
-        else if (e[3] <= e0[ie]) {
-          for (ii = 0; ii < 4; ii++)
-            w1[ib][ie][ii] = 0.25;
-        }
-        else {
-          continue;
-        }
-      }
-    }
-
-    for (ii = 0; ii < 20; ii++)
-      for (ib = 0; ib < nb; ib++)
-        for (ie = 0; ie < ne; ie++)
-          for (jj = 0; jj < 4; jj++)
-            wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (ie = 0; ie < ne; ie++)
-        wght[ik][ib][ie] /= (6.0 * (double)nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-}
diff --git a/src/c/libtetrabz_fermigr.c b/src/c/libtetrabz_fermigr.c
deleted file mode 100644 (file)
index 981b2f7..0000000
+++ /dev/null
@@ -1,512 +0,0 @@
-/*
- Copyright (c) 2014 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include "libtetrabz_common.h"
-#include <stdlib.h>
-
-static void libtetrabz_fermigr3(
-    int ne,
-    double *e0,
-    double de[4],
-    double **w1
-) {
-  int i4, j3, ie, indx[4];
-  double e[4], tsmall[4][3], v;
-
-  for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
-  eig_sort(4, e, indx);
-
-  for (ie = 0; ie < ne; ie++) {
-
-    for (i4 = 0; i4 < 4; i4++) w1[i4][ie] = 0.0;
-
-    if (e[0] < e0[ie] && e0[ie] <= e[1]) {
-
-      libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j3 = 0; j3 < 3; j3++)
-          w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
-    }
-    else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
-
-      libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j3 = 0; j3 < 3; j3++)
-          w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
-
-      libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j3 = 0; j3 < 3; j3++)
-          w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
-    }
-    else if (e[2] < e0[ie] && e0[ie] < e[3]) {
-
-      libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
-      for (i4 = 0; i4 < 4; i4++)
-        for (j3 = 0; j3 < 3; j3++)
-          w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
-    }
-  }
-}
-
-static void libtetrabz_fermigr2(
-    int nb,
-    int ne,
-    double *e0,
-    double *ei1,
-    double **ej1,
-    double ***w1
-) {
-  int ib, i4, j4, ie, indx[4];
-  double e[4], tsmall[4][4], v, de[4], thr, **w2;
-
-  w2 = (double**)malloc(4 * sizeof(double*));
-  w2[0] = (double*)malloc(4 * ne * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    w2[i4] = w2[0] + i4 * ne;
-  }
-
-  thr = 1.0e-8;
-
-  for (ib = 0; ib < nb; ib++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ie = 0; ie < ne; ie++)
-        w1[ib][i4][ie] = 0.0;
-
-    for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
-    eig_sort(4, e, indx);
-
-    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
-
-      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-    }
-    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
-
-      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-
-      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-
-      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-    }
-    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
-
-      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-      
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-
-      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-
-      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_fermigr3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
-      }
-    }
-    else if (e[3] <= 0.0) {
-      for (i4 = 0; i4 < 4; i4++)
-        de[i4] = ej1[ib][i4] - ei1[i4];
-      libtetrabz_fermigr3(ne, e0, de, w2);
-      for (i4 = 0; i4 < 4; i4++)
-        for (ie = 0; ie < ne; ie++)
-          w1[ib][i4][ie] += w2[i4][ie];
-    }
-  }
-
-  free(w2[0]);
-  free(w2);
-}
-
-void fermigr(
-    int ng[3],
-    int nk,
-    int nb,
-    int ne,
-    double bvec[3][3],
-    double **eig1,
-    double **eig2,
-    double *e0,
-    double ****wght
-)
-{
-  /*
-  Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
-  */
-  int it, ik, ib, i20, i4, j4, jb, ie, **ikv, indx[4];
-  double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], **** w1, *** w2, v, tsmall[4][4], thr;
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (it = 0; it < 6 * nk; it++) {
-    ikv[it] = ikv[0] + it * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ej1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  ej1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    ei1[i4] = ei1[0] + i4 * nb;
-    ej1[i4] = ej1[0] + i4 * nb;
-  }
-
-  ej2 = (double**)malloc(nb * sizeof(double*));
-  ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++)
-    ej2[ib] = ej2[0] + ib * 4;
-
-  w1 = (double****)malloc(nb * sizeof(double***));
-  w1[0] = (double***)malloc(nb * 4 * sizeof(double**));
-  w1[0][0] = (double**)malloc(nb * 4 * nb * sizeof(double*));
-  w1[0][0][0] = (double*)malloc(nb * 4 * nb * ne * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
-      for (jb = 0; jb < nb; jb++) {
-        w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
-      }
-    }
-  }
-
-  w2 = (double***)malloc(nb * sizeof(double**));
-  w2[0] = (double**)malloc(nb * 4 * sizeof(double*));
-  w2[0][0] = (double*)malloc(nb * 4 * ne * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w2[ib] = w2[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
-    }
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        for (ie = 0; ie < ne; ie++)
-          wght[ik][ib][jb][ie] = 0.0;
-
-  thr = 1.0e-10;
-
-  for (it = 0; it < 6*nk; it++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ib = 0; ib < nb; ib++) {
-        ei1[i4][ib] = 0.0;
-        ej1[i4][ib] = 0.0;
-      }
-    for (i20 = 0; i20 < 20; i20++) {
-      for (i4 = 0; i4 < 4; i4++) {
-        for (ib = 0; ib < nb; ib++) {
-          ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
-          ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
-        }
-      }
-    }
-    
-    for (ib = 0; ib < nb; ib++) {
-
-      for (i4 = 0; i4 < 4; i4++)
-        for (jb = 0; jb < nb; jb++)
-          for (ie = 0; ie < ne; ie++)
-            w1[ib][i4][jb][ie] = 0.0;
-
-      for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] <= 0.0 && 0.0 < e[1]) {
-
-        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-      }
-      else if (e[1] <= 0.0 && 0.0 < e[2]) {
-
-        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-
-        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-
-        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-      }
-      else if (e[2] <= 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-
-        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-
-        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
-        }
-      }
-      else if (e[3] <= 0.0) {
-        for (i4 = 0; i4 < 4; i4++) {
-          ei2[i4] = ei1[i4][ib];
-          for (jb = 0; jb < nb; jb++)
-            ej2[jb][i4] = ej1[i4][jb];
-        }
-        libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            for (ie = 0; ie < ne; ie++)
-              w1[ib][i4][jb][ie] += w2[jb][i4][ie];
-      }
-      else {
-        continue;
-      }
-    }
-    for (i20 = 0; i20 < 20; i20++)
-      for (ib = 0; ib < nb; ib++)
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            for (ie = 0; ie < ne; ie++)
-              wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][i4] * w1[ib][i4][jb][ie];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        for (ie = 0; ie < ne; ie++)
-          wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(ej1[0]);
-  free(ej1);
-  free(ej2[0]);
-  free(ej2);
-  free(w1[0][0][0]);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-  free(w2[0][0]);
-  free(w2[0]);
-  free(w2);
-}
diff --git a/src/c/libtetrabz_occ.c b/src/c/libtetrabz_occ.c
deleted file mode 100644 (file)
index ca92e9f..0000000
+++ /dev/null
@@ -1,218 +0,0 @@
-//
-// Copyright (c) 2014 Mitsuaki Kawamura
-//
-// Permission is hereby granted, free of charge, to any person obtaining a
-// copy of this software and associated documentation files (the
-// "Software"), to deal in the Software without restriction, including
-// without limitation the rights to use, copy, modify, merge, publish,
-// distribute, sublicense, and/or sell copies of the Software, and to
-// permit persons to whom the Software is furnished to do so, subject to
-// the following conditions:
-// 
-// The above copyright notice and this permission notice shall be included
-// in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-//
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include "libtetrabz_common.h"
-
-void occ(
-  int ng[3],
-  int nk,
-  int nb,
-  double bvec[3][3],
-  double** eig,
-  double** wght
-) {
-  /*
-  Main SUBROUTINE for occupation : Theta(EF - E1)
-  */
-  int it, ik, ib, ii, jj, ** ikv, indx[4];
-  double wlsm[20][4], ** ei1, e[4], ** w1, v, tsmall[4][4];
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (ik = 0; ik < 6 * nk; ik++) {
-    ikv[ik] = ikv[0] + ik * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (ii = 0; ii < 4; ii++) {
-    ei1[ii] = ei1[0] + ii * nb;
-  }
-
-  w1 = (double**)malloc(nb * sizeof(double*));
-  w1[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      wght[ik][ib] = 0.0;
-
-  for (it = 0; it < 6 * nk; it++) {
-
-    for (ii = 0; ii < 4; ii++)
-      for (ib = 0; ib < nb; ib++) 
-        ei1[ii][ib] = 0.0;
-    for (jj = 0; jj < 20; jj++) 
-      for (ii = 0; ii < 4; ii++) 
-        for (ib = 0; ib < nb; ib++) 
-          ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
-
-    for (ib = 0; ib < nb; ib++) {
-
-      for (ii = 0; ii < 4; ii++)
-        w1[ib][ii] = 0.0;
-
-      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] <= 0.0 && 0.0 < e[1]) {
-        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-      }
-      else if (e[1] <= 0.0 && 0.0 < e[2]) {
-
-        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-      }
-      else if (e[2] <= 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-
-        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-        for (ii = 0; ii < 4; ii++)
-          for (jj = 0; jj < 4; jj++)
-            w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
-      }
-      else if (e[3] <= 0.0) {
-        for (ii = 0; ii < 4; ii++)
-          w1[ib][ii] += 0.25;
-      }
-      else {
-        continue;
-      }
-    }
-    for (ii = 0; ii < 20; ii++)
-      for (ib = 0; ib < nb; ib++)
-        for (jj = 0; jj < 4; jj++) {
-          wght[ikv[it][ii]][ib] += wlsm[ii][jj] * w1[ib][jj];
-        }
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      wght[ik][ib] /= (6.0 * (double)nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(w1[0]);
-  free(w1);
-}
-
-void fermieng(
-    int ng[3],
-    int nk,
-    int nb,
-    double bvec[3][3],
-    double **eig,
-    double **wght,
-  double nelec,
-  int *iteration,
-  double *ef
-) {
-  /*
-  Calculate Fermi energy
-  */
-  int maxiter, ik, ib, iter;
-  double eps, elw, eup, **eig2, sumkmid;
-
-  eig2 = (double**)malloc(nk * sizeof(double*));
-  eig2[0] = (double*)malloc(nk * nb * sizeof(double));
-  for (ik = 0; ik < nk; ik++) {
-    eig2[ik] = eig2[0] + ik * nb;
-  }
-
-  maxiter = 300;
-  eps = 1.0e-10;
-
-  elw = eig[0][0];
-  eup = eig[0][0];
-  for (ik = 0; ik < nk; ik++) {
-    for (ib = 0; ib < nb; ib++) {
-      if (elw > eig[ik][ib]) elw = eig[ik][ib];
-      if (eup < eig[ik][ib]) eup = eig[ik][ib];
-    }
-  }
-  /*
-  Bisection method
-  */
-  for (iter = 0; iter < maxiter; iter++) {
-
-    *ef = (eup + elw) * 0.5;
-    /*
-    Calc. # of electrons
-    */
-    for (ik = 0; ik < nk; ik++)
-      for (ib = 0; ib < nb; ib++)
-        eig2[ik][ib] = eig[ik][ib] - *ef;
-    occ(ng, nk, nb, bvec, eig2, wght);
-
-    sumkmid = 0.0;
-    for (ik = 0; ik < nk; ik++) {
-      for (ib = 0; ib < nb; ib++) {
-        sumkmid += wght[ik][ib];
-      }
-    }
-    /*
-    convergence check
-    */
-    if (fabs(sumkmid - nelec) < eps) {
-      *iteration = iter;
-      free(eig2);
-      return;
-    }
-    else if (sumkmid < nelec)
-      elw = *ef;
-    else
-      eup = *ef;
-  }
-  free(eig2);
-}
diff --git a/src/c/libtetrabz_polcmplx.c b/src/c/libtetrabz_polcmplx.c
deleted file mode 100644 (file)
index 88b53fa..0000000
+++ /dev/null
@@ -1,906 +0,0 @@
-/*
- Copyright (c) 2014 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
-
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include "libtetrabz_common.h"
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-
-static void libtetrabz_polcmplx_1234(
-  double g1,
-  double g2,
-  double g3,
-  double g4,
-  double* w
-) {
-  /*
-  1, Different each other
-  */
-  double w2, w3, w4;
-  /*
-  Real
-  */
-  w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
-    (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
-    (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
-    (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
-  w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
-  w4 = w4 / (g4 - g1);
-  w2 = (w2 - w3) / (g2 - g3);
-  w4 = (w4 - w3) / (g4 - g3);
-  w[0] = (w4 - w2) / (2.0 * (g4 - g2));
-  /*
-  Imaginal
-  */
-  w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
-    (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = 4.0 * g2 - w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
-    (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = 4.0 * g3 - w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
-    (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
-  w4 = 4.0 * g4 - w4 / (g4 - g1);
-  w4 = w4 / (g4 - g1);
-  w2 = (w2 - w3) / (g2 - g3);
-  w4 = (w4 - w3) / (g4 - g3);
-  w[1] = (w4 - w2) / (2.0 * (g4 - g2));
-}
-
-
-static void libtetrabz_polcmplx_1231(
-  double g1,
-  double g2,
-  double g3,
-  double* w
-) {
-  /*
-  2, g4 = g1
-  */
-  double w2, w3;
-  /*
-  Real
-  */
-  w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
-    + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
-  w2 = -g1 + w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
-    + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
-  w3 = -g1 + w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w[0] = (w3 - w2) / (2.0 * (g3 - g2));
-  /*
-  Imaginal
-  */
-  w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
-    (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = 4.0 * g2 - w2 / (g2 - g1);
-  w2 = 1 + w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
-    (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = 4.0 * g3 - w3 / (g3 - g1);
-  w3 = 1 + w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w[1] = (w3 - w2) / (2.0 * (g3 - g2));
-}
-
-
-static void libtetrabz_polcmplx_1233(
-  double g1,
-  double g2,
-  double g3,
-  double* w
-) {
-  /*
-  3, g4 = g3
-  */
-  double w2, w3;
-  /*
-  Real
-  */
-  w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
-    g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
-    g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w2 = (w3 - w2) / (g3 - g2);
-  w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
-    + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
-  w3 = 4.0 * g1 + w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w[0] = (w3 - w2) / (2.0 * (g3 - g2));
-  /*
-  Imaginal
-  */
-  w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
-    (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w2 = 4.0 * g2 - w2 / (g2 - g1);
-  w2 = w2 / (g2 - g1);
-  w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
-    (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = 4.0 * g3 - w3 / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w2 = (w3 - w2) / (g3 - g2);
-  w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3 * g3) * (atan(g3) - atan(g1))
-    + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
-  w3 = w3 / (g3 - g1) - 4.0 * g1;
-  w3 = w3 / (g3 - g1) - 2.0;
-  w3 = (2.0 * w3) / (g3 - g1);
-  w[1] = (w3 - w2) / (2.0 * (g3 - g2));
-}
-
-static void libtetrabz_polcmplx_1221(
-  double g1,
-  double g2,
-  double* w
-) {
-  /*
-  4, g4 = g1 and g3 = g2
-  */
-  /*
-  Real
-  */
-  w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
-    + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
-  w[0] = 3.0 * g1 + w[0] / (g2 - g1);
-  w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
-  w[0] = w[0] / (2.0 * (g2 - g1));
-  /*
-  Imaginal
-  */
-  w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
-    + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
-  w[1] = -4.0 * g1 + w[1] / (g2 - g1);
-  w[1] = -3.0 + w[1] / (g2 - g1);
-  w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
-}
-
-
-static void libtetrabz_polcmplx_1222(
-  double g1,
-  double g2,
-  double* w
-) {
-  /*
-  5, g4 = g3 = g2
-  */
-  /*
-  Real
-  */
-  w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
-    + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
-  w[0] = g1 - w[0] / (g2 - g1);
-  w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
-  w[0] = w[0] / (2.0 * (g2 - g1));
-  /*
-  Imaginal
-  */
-  w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
-    + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w[1] = 4.0 * g1 + w[1] / (g2 - g1);
-  w[1] = 1.0 + w[1] / (g2 - g1);
-  w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
-}
-
-
-static void libtetrabz_polcmplx_1211(
-  double g1,
-  double g2,
-  double* w
-) {
-  /*
-  6, g4 = g3 = g1
-  */
-  /*
-  Real
-  */
-  w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
-    + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
-  w[0] = -5.0 * g1 + w[0] / (g2 - g1);
-  w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
-  w[0] = w[0] / (6.0 * (g2 - g1));
-  /*
-  Imaginal
-  */
-  w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
-    + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
-  w[1] = 4.0 * g2 + w[1] / (g2 - g1);
-  w[1] = 1.0 + w[1] / (g2 - g1);
-  w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
-}
-
-static void libtetrabz_polcmplx3(
-  int ne,
-  double** e0,
-  double de[4],
-  double*** w1
-) {
-  /*
-  Tetrahedron method for delta(om - ep + e)
-  */
-  int i4, ir, indx[4], ie;
-  double e[4], x[4], thr, w2[4][2], denom;
-
-  for (i4 = 0; i4 < 3; i4++) e[i4] = de[i4];
-  eig_sort(4, e, indx);
-
-  for (ie = 0; ie < ne; ie++) {
-    /*
-    I am not sure which one is better.
-    The former is more stable.
-    The latter is more accurate ?
-    */
-    for (i4 = 0; i4 < 4; i4++) {
-      denom = (de[i4] + e0[ie][0]) * (de[i4] + e0[ie][0]) + e0[ie][1] * e0[ie][1];
-      w1[i4][ie][0] = 0.25 * (de[i4] + e0[ie][0]) / denom;
-      w1[i4][ie][1] = 0.25 * (-e0[ie][1]) / denom;
-    }
-    continue;
-
-    for (i4 = 0; i4 < 4; i4++)
-      x[i4] = (e[i4] + e0[ie][0]) / e0[ie][1];
-    /* thr = maxval(de(1:4)) * 1d-3*/
-    thr = fabs(x[0]);
-    for (i4 = 0; i4 < 4; i4++)
-      if (thr < fabs(x[i4])) thr = fabs(x[i4]);
-    thr = fmax(1.0e-3, thr * 1.0e-2);
-
-    if (fabs(x[3] - x[2]) < thr) {
-      if (fabs(x[3] - x[1]) < thr) {
-        if (fabs(x[3] - x[0]) < thr) {
-          /*
-          e[3] = e[2] = e[1] = e[0]
-          */
-          w2[3][0] = 0.25 * x[3] / (1.0 + x[3] * x[3]);
-          w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
-          for (ir = 0; ir < 2; ir++) {
-            w2[2][ir] = w2[3][ir];
-            w2[1][ir] = w2[3][ir];
-            w2[0][ir] = w2[3][ir];
-          }
-        }
-        else {
-          /*
-          e[3] = e[2] = e[1]
-          */
-          libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
-          for (ir = 0; ir < 2; ir++) {
-            w2[2][ir] = w2[3][ir];
-            w2[1][ir] = w2[3][ir];
-          }
-          libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
-          /*
-          # IF(ANY(w2(1:2,1:4) < 0.0)):
-          #   WRITE(*,*) ie
-          #   WRITE(*,'(100e15.5)') x[0:4]
-          #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
-          #   STOP "weighting 4=3=2"
-          */
-        }
-      }
-      else if (fabs(x[1] - x[0]) < thr) {
-        /*
-        e[3] = e[2], e[1] = e[0]
-        */
-        libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
-        for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
-        libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
-        for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
-        /*
-         IF(ANY(w2(1:2,1:4) < 0.0)):
-           WRITE(*,*) ie
-           WRITE(*,'(100e15.5)') x[0:4]
-           WRITE(*,'(2e15.5)') w2(1:2,1:4)
-           STOP "weighting 4=3 2=1"
-        */
-      }
-      else {
-        /*
-        e[3] = e[2]
-        */
-        libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
-        for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
-        libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
-        libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
-        /*
-         IF(ANY(w2(1:2,1:4) < 0.0)):
-           WRITE(*,*) ie
-           WRITE(*,'(100e15.5)') x[0:4]
-           WRITE(*,'(2e15.5)') w2(1:2,1:4)
-           STOP "weighting 4=3"
-        */
-      }
-    }
-    else if (fabs(x[2] - x[1]) < thr) {
-      if (fabs(x[2] - x[0]) < thr) {
-        /*
-        e[2] = e[1] = e[0]
-        */
-        libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
-        libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
-        for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
-        for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[2][ir];
-        /*
-        IF(ANY(w2(1:2,1:4) < 0.0)):
-          WRITE(*,*) ie
-          WRITE(*,'(100e15.5)') x[0:4]
-          WRITE(*,'(2e15.5)') w2(1:2,1:4)
-          STOP "weighting 3=2=1"
-        */
-      }
-      else {
-        /*
-        e[2] = e[1]
-        */
-        libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
-        libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
-        for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
-        libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
-        /*
-        IF(ANY(w2(1:2,1:4) < 0.0)):
-          WRITE(*,*) ie
-          WRITE(*,'(100e15.5)') x[0:4]
-          WRITE(*,'(2e15.5)') w2(1:2,1:4)
-          STOP "weighting 3=2"
-        */
-      }
-    }
-    else if (fabs(x[1] - x[0]) < thr) {
-      /*
-      e[1] = e[0]
-      */
-      libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
-      libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
-      libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
-      for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
-      /*
-      IF(ANY(w2(1:2,1:4) < 0.0)):
-        WRITE(*,*) ie
-        WRITE(*,'(100e15.5)') x[0:4]
-        WRITE(*,'(2e15.5)') w2(1:2,1:4)
-        STOP "weighting 2=1"
-      */
-    }
-    else {
-      /*
-      Different each other.
-      */
-      libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2], w2[3]);
-      libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3], w2[2]);
-      libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3], w2[1]);
-      libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3], w2[0]);
-      /*
-      IF(ANY(w2(1:2,1:4) < 0.0)):
-        WRITE(*,*) ie
-        WRITE(*,'(100e15.5)') x[0:4]
-        WRITE(*,'(2e15.5)') w2(1:2,1:4)
-        STOP "weighting"
-      */
-    }
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[indx[i4]][ie][0] = w2[i4][0] / e0[ie][1];
-      w1[indx[i4]][ie][1] = w2[i4][1] / (-e0[ie][1]);
-    }
-  }
-}
-//
-// Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
-//  for 0<x<1, 0<y<1-x, 0<z<1-x-y
-
-static void libtetrabz_polcmplx2(
-  int nb,
-  int ne,
-  double* *e0,
-  double* ei1,
-  double** ej1,
-  double**** w1
-) {
-  /*
-  Tetrahedron method for theta( - E2)
-  */
-  int ib, i4, j4, ir, ie, indx[4];
-  double e[4], tsmall[4][4], v, de[4], thr, *** w2;
-  
-  w2 = (double***)malloc(4 * sizeof(double**));
-  w2[0] = (double**)malloc(4 * ne * sizeof(double*));
-  w2[0][0] = (double*)malloc(4 * ne * 2 * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    w2[i4] = w2[0] + i4 * ne;
-    for (ie = 0; ie < ne; ie++) {
-      w2[i4][ie] = w2[0][0] + i4 * ne * 2 + ie * 2;
-    }
-  }
-
-  thr = 1.0e-8;
-
-  for (ib = 0; ib < nb; ib++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ie = 0; ie < ne; ie++)
-        for (ir = 0; ir < 2; ir++)
-          w1[ib][i4][ie][ir] = 0.0;
-
-    for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
-    eig_sort(4, e, indx);
-
-    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
-
-      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-    }
-    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
-
-      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-
-      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-
-      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-    }
-    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
-
-      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-
-      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-
-      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polcmplx3(ne, e0, de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
-      }
-    }
-    else if (e[3] <= 0.0) {
-      for (i4 = 0; i4 < 4; i4++)
-        de[i4] = ej1[ib][i4] - ei1[i4];
-      libtetrabz_polcmplx3(ne, e0, de, w2);
-      for (i4 = 0; i4 < 4; i4++)
-        for (ie = 0; ie < ne; ie++)
-          for (ir = 0; ir < 2; ir++)
-            w1[ib][i4][ie][ir] += w2[i4][ie][ir];
-    }
-  }
-  free(w2[0][0]);
-  free(w2[0]);
-  free(w2);
-}
-
-
-void polcmplx(
-    int ng[3],
-    int nk,
-    int nb,
-    int ne,
-    double bvec[3][3],
-    double **eig1,
-    double **eig2,
-    double **e0,
-    double *****wght) {
-  /*
-  Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
-  */
-  int it, ik, ie, ib, i4, j4, ir, jb, **ikv, indx[4], i20;
-  double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], ***** w1, **** w2, v, tsmall[4][4], thr;
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (it = 0; it < 6 * nk; it++) {
-    ikv[it] = ikv[0] + it * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ej1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  ej1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    ei1[i4] = ei1[0] + i4 * nb;
-    ej1[i4] = ej1[0] + i4 * nb;
-  }
-
-  ej2 = (double**)malloc(nb * sizeof(double*));
-  ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++)
-    ej2[ib] = ej2[0] + ib * 4;
-
-  w1 = (double*****)malloc(nb * sizeof(double****));
-  w1[0] = (double****)malloc(nb * 4 * sizeof(double***));
-  w1[0][0] = (double***)malloc(nb * 4 * nb * sizeof(double**));
-  w1[0][0][0] = (double**)malloc(nb * 4 * nb * ne * sizeof(double*));
-  w1[0][0][0][0] = (double*)malloc(nb * 4 * nb * ne * 2 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
-      for (jb = 0; jb < nb; jb++) {
-        w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
-        for (ie = 0; ie < ne; ie++) {
-          w1[ib][i4][jb][ie] = w1[0][0][0][0] + ib * 4 * nb * ne * 2 + i4 * nb * ne * 2 + jb * ne * 2 + ie * 2;
-        }
-      }
-    }
-  }
-
-  w2 = (double****)malloc(nb * sizeof(double***));
-  w2[0] = (double***)malloc(nb * 4 * sizeof(double**));
-  w2[0][0] = (double**)malloc(nb * 4 * ne * sizeof(double*));
-  w2[0][0][0] = (double*)malloc(nb * 4 * ne * 2 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w2[ib] = w2[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
-      for (ie = 0; ie < ne; ie++) {
-        w2[ib][i4][ie] = w2[0][0][0] + ib * 4 * ne * 2 + i4 * ne * 2 + ie * 2;
-      }
-    }
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        for (ie = 0; ie < ne; ie++)
-          for (ir = 0; ir < 2; ir++)
-            wght[ik][ib][jb][ie][ir] = 0.0;
-
-  thr = 1.0e-8;
-
-  for (it = 0; it < 6 * nk; it++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ib = 0; ib < nb; ib++) {
-        ei1[i4][ib] = 0.0;
-        ej1[i4][ib] = 0.0;
-      }
-    for (i20 = 0; i20 < 20; i20++) {
-      for (i4 = 0; i4 < 4; i4++) {
-        for (ib = 0; ib < nb; ib++) {
-          ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
-          ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
-        }
-      }
-    }
-
-    for (ib = 0; ib < nb; ib++) {
-
-      for (i4 = 0; i4 < 4; i4++)
-        for (jb = 0; jb < nb; jb++)
-          for (ie = 0; ie < ne; ie++)
-            for (ir = 0; ir < 2; ir++)
-              w1[ib][i4][jb][ie][ir] = 0.0;
-
-      for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] <= 0.0 && 0.0 < e[1]) {
-
-        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-      }
-      else if (e[1] <= 0.0 && 0.0 < e[2]) {
-
-        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-
-        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-
-        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-      }
-      else if (e[2] <= 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-
-        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-
-        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                for (ie = 0; ie < ne; ie++)
-                  for (ir = 0; ir < 2; ir++)
-                    w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
-        }
-      }
-      else if (e[3] <= 0.0) {
-        for (i4 = 0; i4 < 4; i4++) {
-          ei2[i4] = ei1[i4][ib];
-          for (jb = 0; jb < nb; jb++)
-            ej2[jb][i4] = ej1[i4][jb];
-        }
-        libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                w1[ib][i4][jb][ie][ir] += w2[jb][i4][ie][ir];
-      }
-      else {
-        continue;
-      }
-    }
-    for (i20 = 0; i20 < 20; i20++)
-      for (ib = 0; ib < nb; ib++)
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            for (ie = 0; ie < ne; ie++)
-              for (ir = 0; ir < 2; ir++)
-                wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][i4] * w1[ib][i4][jb][ie][ir];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        for (ie = 0; ie < ne; ie++)
-          for (i4 = 0; i4 < 2; i4++)
-            wght[ik][ib][jb][ie][i4] /= (6.0 * (double) nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(ej1[0]);
-  free(ej1);
-  free(ej2[0]);
-  free(ej2);
-  free(w1[0][0][0][0]);
-  free(w1[0][0][0]);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-  free(w2[0][0][0]);
-  free(w2[0][0]);
-  free(w2[0]);
-  free(w2);
-}
diff --git a/src/c/libtetrabz_polstat.c b/src/c/libtetrabz_polstat.c
deleted file mode 100644 (file)
index a56b327..0000000
+++ /dev/null
@@ -1,705 +0,0 @@
-/*
- Copyright (c) 2014 Mitsuaki Kawamura
-
- Permission is hereby granted, free of charge, to any person obtaining a
- copy of this software and associated documentation files (the
- "Software"), to deal in the Software without restriction, including
- without limitation the rights to use, copy, modify, merge, publish,
- distribute, sublicense, and/or sell copies of the Software, and to
- permit persons to whom the Software is furnished to do so, subject to
- the following conditions:
- The above copyright notice and this permission notice shall be included
- in all copies or substantial portions of the Software.
-
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#include "libtetrabz_common.h"
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-
-/*
- Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
- for 0<x<1, 0<y<1-x, 0<z<1-x-y
-*/
-static double libtetrabz_polstat_1234(
-    double g1,
-    double g2,
-    double g3,
-    double g4,
-    double lng1,
-    double lng2,
-    double lng3,
-    double lng4)
-{
-  /*
-  1, Different each other
-  */
-  double w2, w3, w4, w;
-  w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 / (g2 - g1);
-  w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 / (g3 - g1);
-  w4 = ((lng4 - lng1) / (g4 - g1) * g4 - 1.0) * g4 / (g4 - g1);
-  w2 = ((w2 - w3) * g2) / (g2 - g3);
-  w4 = ((w4 - w3) * g4) / (g4 - g3);
-  w = (w4 - w2) / (g4 - g2);
-  return w;
-}
-
-
-static double libtetrabz_polstat_1231(
-    double g1,
-    double g2,
-    double g3,
-    double lng1,
-    double lng2,
-    double lng3
-) {
-  /*
-  2, g4 = g1
-  */
-  double w2, w3, w;
-  w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 * g2 / (g2 - g1) - g1 / 2.0;
-  w2 = w2 / (g2 - g1);
-  w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 * g3 / (g3 - g1) - g1 / 2.0;
-  w3 = w3 / (g3 - g1);
-  w = (w3 - w2) / (g3 - g2);
-  return w;
-}
-
-static double libtetrabz_polstat_1233(
-    double g1,
-    double g2,
-    double g3,
-    double lng1,
-    double lng2,
-    double lng3
-) {
-  /*
-  # 3, g4 = g3
-  */
-  double w2, w3, w;
-  w2 = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
-  w2 = (g2 * w2) / (g2 - g1);
-  w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
-  w3 = (g3 * w3) / (g3 - g1);
-  w2 = (w3 - w2) / (g3 - g2);
-  w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
-  w3 = 1.0 - (2.0 * w3 * g1) / (g3 - g1);
-  w3 = w3 / (g3 - g1);
-  w = (g3 * w3 - g2 * w2) / (g3 - g2);
-  return w;
-}
-
-static double libtetrabz_polstat_1221(
-    double g1,
-    double g2,
-    double lng1,
-    double lng2
-) {
-  /*
-  4, g4 = g1 and g3 = g2
-  */
-  double w;
-  w = 1.0 - (lng2 - lng1) / (g2 - g1) * g1;
-  w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
-  w = -1.0 + (3.0 * g2 * w) / (g2 - g1);
-  w = w / (2.0 * (g2 - g1));
-  return w;
-}
-
-static double libtetrabz_polstat_1222(
-    double g1,
-    double g2,
-    double lng1,
-    double lng2
-) {
-  /*
-  5, g4 = g3 = g2
-  */
-  double w;
-  w = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
-  w = (2.0 * g1 * w) / (g2 - g1) - 1.0;
-  w = (3.0 * g1 * w) / (g2 - g1) + 1.0;
-  w = w / (2.0 * (g2 - g1));
-  return w;
-}
-
-static double libtetrabz_polstat_1211(
-    double g1,
-    double g2,
-    double lng1,
-    double lng2
-) {
-  /*
-  6, g4 = g3 = g1
-  */
-  double w;
-  w = -1.0 + (lng2 - lng1) / (g2 - g1) * g2;
-  w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
-  w = -1.0 + (3.0 * g2 * w) / (2.0 * (g2 - g1));
-  w = w / (3.0 * (g2 - g1));
-  return w;
-}
-
-static void libtetrabz_polstat3(
-    double de[4],
-    double w1[4]
-)
-{
-  /*
-  Tetrahedron method for delta(om - ep + e)
-  */
-  int i4, indx[4];
-  double thr, thr2, e[4], ln[4];
-
-  for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
-  eig_sort(4, e, indx);
-
-  thr = e[3] * 1.0e-3;
-  thr2 = 1.0e-8;
-
-  for(i4 =0; i4 <4; i4++){
-    if (e[i4] < thr2) {
-      if (i4 == 3) {
-        printf("  Nesting # \n");
-      }
-      ln[i4] = 0.0;
-      e[i4] = 0.0;
-    }
-    else{
-      ln[i4] = log(e[i4]);
-    }
-  }
-
-  if (fabs(e[3] - e[2]) < thr) {
-    if (fabs(e[3] - e[1]) < thr) {
-      if (fabs(e[3] - e[0]) < thr) {
-        /*
-        e[3] = e[2] = e[1] = e[0]
-        */
-        w1[indx[3]] = 0.25 / e[3];
-        w1[indx[2]] = w1[indx[3]];
-        w1[indx[1]] = w1[indx[3]];
-        w1[indx[0]] = w1[indx[3]];
-      }
-      else {
-        /*
-        e[3] = e[2] = e[1]
-        */
-        w1[indx[3]] = libtetrabz_polstat_1211(e[3], e[0], ln[3], ln[0]);
-        w1[indx[2]] = w1[indx[3]];
-        w1[indx[1]] = w1[indx[3]];
-        w1[indx[0]] = libtetrabz_polstat_1222(e[0], e[3], ln[0], ln[3]);
-        if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-          printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-          printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-          printf("weighting 4=3=2\n");
-        }
-      }
-    }
-    else if (fabs(e[1] - e[0]) < thr) {
-      /*
-      e[3] = e[2], e[1] = e[0]
-      */
-      w1[indx[3]] = libtetrabz_polstat_1221(e[3], e[1], ln[3], ln[1]);
-      w1[indx[2]] = w1[indx[3]];
-      w1[indx[1]] = libtetrabz_polstat_1221(e[1], e[3], ln[1], ln[3]);
-      w1[indx[0]] = w1[indx[1]];
-
-      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-        printf("weighting 4=3 2=1\n");
-      }
-    }
-    else {
-      /*
-      e[3] = e[2]
-      */
-      w1[indx[3]] = libtetrabz_polstat_1231(e[3], e[0], e[1], ln[3], ln[0], ln[1]);
-      w1[indx[2]] = w1[indx[3]];
-      w1[indx[1]] = libtetrabz_polstat_1233(e[1], e[0], e[3], ln[1], ln[0], ln[3]);
-      w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[1], e[3], ln[0], ln[1], ln[3]);
-
-      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-        printf("weighting 4=3\n");
-      }
-    }
-  }  
-  else if (fabs(e[2] - e[1]) < thr) {
-    if (fabs(e[2] - e[0]) < thr) {
-      /*
-      e[2] = e[1] = e[0]
-      */
-      w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
-      w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
-      w1[indx[1]] = w1[indx[2]];
-      w1[indx[0]] = w1[indx[2]];
-
-      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-        printf("weighting 3=2=1\n");
-      }
-    }
-    else {
-      /*
-      e[2] = e[1]
-      */
-      w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
-      w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
-      w1[indx[1]] = w1[indx[2]];
-      w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
-
-      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-        printf("weighting 3=2\n");
-      }
-    }
-  }  
-  else if (fabs(e[1] - e[0]) < thr) {
-    /*
-    e[1] = e[0]
-    */
-    w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
-    w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
-    w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
-    w1[indx[0]] = w1[indx[1]];
-
-    if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-      printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-      printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-      printf("weighting 2=1\n");
-    }
-  }  
-  else {
-    /*
-    Different each other.
-    */
-    w1[indx[3]] = libtetrabz_polstat_1234(e[3], e[0], e[1], e[2], ln[3], ln[0], ln[1], ln[2]);
-    w1[indx[2]] = libtetrabz_polstat_1234(e[2], e[0], e[1], e[3], ln[2], ln[0], ln[1], ln[3]);
-    w1[indx[1]] = libtetrabz_polstat_1234(e[1], e[0], e[2], e[3], ln[1], ln[0], ln[2], ln[3]);
-    w1[indx[0]] = libtetrabz_polstat_1234(e[0], e[1], e[2], e[3], ln[0], ln[1], ln[2], ln[3]);
-
-    if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
-      printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
-      printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
-      printf("weighting\n");
-    }
-  }
-}
-
-static void libtetrabz_polstat2(
-    int nb,
-    double *ei1,
-    double **ej1,
-    double **w1
-) {
-  /*
-  Tetrahedron method for theta( - E2)
-  :param ei1:
-  :param ej1:
-  :return:
-  */
-  int i4, j4, ib, indx[4];
-  double v, thr, e[4], de[4], w2[4], tsmall[4][4];
-
-  thr = 1.0e-8;
-
-  for (ib = 0; ib < nb; ib++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      w1[ib][i4] = 0.0;
-
-    for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
-    eig_sort(4, e, indx);
-
-    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
-
-      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-    }
-    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
-
-      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-
-      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-
-      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-    }
-    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
-
-      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-      
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-
-      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-
-      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-      if (v > thr) {
-        for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
-        libtetrabz_polstat3(de, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (j4 = 0; j4 < 4; j4++)
-            w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
-      }
-    }
-    else if (e[3] <= 0.0) {
-      for (i4 = 0; i4 < 4; i4++)
-        de[i4] = ej1[ib][i4] - ei1[i4];
-      libtetrabz_polstat3(de, w2);
-      for (i4 = 0; i4 < 4; i4++) 
-        w1[ib][i4] += w2[i4];
-    }
-  }
-}
-
-void polstat(
-    int ng[3],
-    int nk,
-    int nb,
-    double bvec[3][3],
-    double **eig1,
-    double **eig2,
-    double ***wght
-)
-{
-    /*
-    Compute Static polarization function
-    */
-  int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4];
-  double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr;
-
-  thr = 1.0e-10;
-
-  ikv = (int**)malloc(6 * nk * sizeof(int*));
-  ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
-  for (it = 0; it < 6 * nk; it++) {
-    ikv[it] = ikv[0] + it * 20;
-  }
-
-  ei1 = (double**)malloc(4 * sizeof(double*));
-  ej1 = (double**)malloc(4 * sizeof(double*));
-  ei1[0] = (double*)malloc(4 * nb * sizeof(double));
-  ej1[0] = (double*)malloc(4 * nb * sizeof(double));
-  for (i4 = 0; i4 < 4; i4++) {
-    ei1[i4] = ei1[0] + i4 * nb;
-    ej1[i4] = ej1[0] + i4 * nb;
-  }
-
-  ej2 = (double**)malloc(nb * sizeof(double*));
-  ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++)
-    ej2[ib] = ej2[0] + ib * 4;
-
-  w1 = (double***)malloc(nb * sizeof(double**));
-  w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
-  w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w1[ib] = w1[0] + ib * 4;
-    for (i4 = 0; i4 < 4; i4++) {
-      w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
-    }
-  }
-
-  w2 = (double**)malloc(nb * sizeof(double*));
-  w2[0] = (double*)malloc(nb * 4 * sizeof(double));
-  for (ib = 0; ib < nb; ib++) {
-    w2[ib] = w2[0] + ib * 4;
-  }
-
-  libtetrabz_initialize(ng, bvec, wlsm, ikv);
-
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] = 0.0;
-
-  for(it = 0; it < 6*nk; it++) {
-
-    for (i4 = 0; i4 < 4; i4++)
-      for (ib = 0; ib < nb; ib++) {
-        ei1[i4][ib] = 0.0;
-        ej1[i4][ib] = 0.0;
-      }
-    for (i20 = 0; i20 < 20; i20++) {
-      for (i4 = 0; i4 < 4; i4++) {
-        for (ib = 0; ib < nb; ib++) {
-          ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
-          ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
-        }
-      }
-    }
-    
-    for (ib = 0; ib < nb; ib++) {
-
-      for (i4 = 0; i4 < 4; i4++)
-        for (jb = 0; jb < nb; jb++)
-          w1[ib][i4][jb] = 0.0;
-
-      for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
-      eig_sort(4, e, indx);
-
-      if (e[0] <= 0.0 && 0.0 < e[1]) {
-
-        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[1] <= 0.0 && 0.0 < e[2]) {
-
-        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[2] <= 0.0 && 0.0 < e[3]) {
-
-        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-
-        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
-
-        if (v > thr) {
-          for (j4 = 0; j4 < 4; j4++) {
-            ei2[j4] = 0.0;
-            for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
-          }
-          for (i4 = 0; i4 < 4; i4++)
-            for (j4 = 0; j4 < 4; j4++) {
-              ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
-              for (jb = 0; jb < nb; jb++)
-                ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
-            }
-          libtetrabz_polstat2(nb, ei2, ej2, w2);
-          for (i4 = 0; i4 < 4; i4++)
-            for (jb = 0; jb < nb; jb++)
-              for (j4 = 0; j4 < 4; j4++)
-                w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
-        }
-      }
-      else if (e[3] <= 0.0) {
-        for (i4 = 0; i4 < 4; i4++) {
-          ei2[i4] = ei1[i4][ib];
-          for (jb = 0; jb < nb; jb++)
-            ej2[jb][i4] = ej1[i4][jb];
-        }
-        libtetrabz_polstat2(nb, ei2, ej2, w2);
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            w1[ib][i4][jb] += w2[jb][i4];
-      }
-      else {
-        continue;
-      }
-    }
-    for (i20 = 0; i20 < 20; i20++)
-      for (ib = 0; ib < nb; ib++)
-        for (i4 = 0; i4 < 4; i4++)
-          for (jb = 0; jb < nb; jb++)
-            wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
-  }
-  for (ik = 0; ik < nk; ik++)
-    for (ib = 0; ib < nb; ib++)
-      for (jb = 0; jb < nb; jb++)
-        wght[ik][ib][jb] /= (6.0 * (double) nk);
-
-  free(ikv[0]);
-  free(ikv);
-  free(ei1[0]);
-  free(ei1);
-  free(ej1[0]);
-  free(ej1);
-  free(ej2[0]);
-  free(ej2);
-  free(w1[0][0]);
-  free(w1[0]);
-  free(w1);
-  free(w2[0]);
-  free(w2);
-}
diff --git a/src/c/test.c b/src/c/test.c
deleted file mode 100644 (file)
index 1b13469..0000000
+++ /dev/null
@@ -1,682 +0,0 @@
-/*\r
-# Copyright (c) 2014 Mitsuaki Kawamura\r
-#\r
-# Permission is hereby granted, free of charge, to any person obtaining a\r
-# copy of this software and associated documentation files (the\r
-# "Software"), to deal in the Software without restriction, including\r
-# without limitation the rights to use, copy, modify, merge, publish,\r
-# distribute, sublicense, and/or sell copies of the Software, and to\r
-# permit persons to whom the Software is furnished to do so, subject to\r
-# the following conditions:\r
-#\r
-# The above copyright notice and this permission notice shall be included\r
-# in all copies or substantial portions of the Software.\r
-#\r
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS\r
-# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
-# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.\r
-# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY\r
-# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,\r
-# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE\r
-# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
-*/\r
-#include <stdio.h>\r
-#include <stdlib.h>\r
-#include <math.h>\r
-#include "libtetrabz.h"\r
-\r
-const double pi = 3.1415926535897932384626433832795;\r
-\r
-\r
-static void test_occ(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig,\r
-  double* mat\r
-) {\r
-  int ib, ik;\r
-  double val[2], ** wght, **eig2;\r
-\r
-  wght = (double**)malloc(nk * sizeof(double*));\r
-  eig2 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig2[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig2[ik] = eig2[0] + ik * nb;\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      eig2[ik][ib] = eig[ik][ib] - 0.5;\r
-\r
-  occ(ng, nk, nb, bvec, eig2, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++) val[ib] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      val[ib] += wght[ik][ib] * mat[ik];\r
-  printf("# libtetrabz_occ\n");\r
-  printf("     % 15.5e % 15.5e\n", 4.0 * pi / 5.0, val[0] * vbz);\r
-  printf("     % 15.5e % 15.5e\n", pi / (5.0 * sqrt(2.0)), val[1] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig2[0]);\r
-  free(eig2);\r
-}\r
-\r
-static void test_fermieng(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig,\r
-  double* mat\r
-) {\r
-  int ib, ik, iteration;\r
-  double val[2], ** wght, nelec, ef;\r
-\r
-  nelec = (4.0 * pi / 3.0 + sqrt(2.0) * pi / 3.0) / vbz;\r
-\r
-  wght = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-  }\r
-\r
-  fermieng(ng, nk, nb, bvec, eig, wght, nelec, &iteration, &ef);\r
-  printf("debug %d\n", iteration);\r
-\r
-  for (ib = 0; ib < nb; ib++) val[ib] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      val[ib] += wght[ik][ib] * mat[ik];\r
-\r
-  printf("# libtetrabz_fermieng\n");\r
-  printf("     %15.5e %15.5e\n", 0.5, ef);\r
-  printf("     %15.5e %15.5e\n", 4.0 * pi / 5.0, val[0] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi / (5.0 * sqrt(2.0)), val[1] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0]);\r
-  free(wght);\r
-}\r
-\r
-\r
-static void test_dos(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig,\r
-  double* mat\r
-) {\r
-  int ib, ik, ne, ie;\r
-  double *** wght, e0[5], val0[2][5], val[2][5];\r
-\r
-  ne = 5;\r
-  wght = (double***)malloc(nk * sizeof(double**));\r
-  wght[0] = (double**)malloc(nk * nb * sizeof(double*));\r
-  wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;\r
-    }\r
-  }\r
-\r
-  for (ie = 0; ie < ne; ie++) {\r
-    e0[ie] = 0.2 * (double)(ie + 1);\r
-    val0[0][ie] = 4.0 * pi * e0[ie] * e0[ie] * e0[ie];\r
-    if (e0[ie] > 1.0 / sqrt(2.0))\r
-      val0[1][ie] = sqrt(2.0) * pi * pow(sqrt(-1.0 + 2.0 * e0[ie] * e0[ie]), 3);\r
-    else\r
-      val0[1][ie] = 0.0;\r
-    e0[ie] = e0[ie] * e0[ie] * 0.5;\r
-  }\r
-\r
-  dos(ng, nk, nb, ne, bvec, eig, e0, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (ie = 0; ie < ne; ie++)\r
-      val[ib][ie] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (ie = 0; ie < ne; ie++)\r
-        val[ib][ie] += wght[ik][ib][ie] * mat[ik];\r
-\r
-  printf("# libtetrabz_dos\n");\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (ie = 0; ie < ne; ie++)\r
-      printf("     %15.5e %15.5e\n", val0[ib][ie], val[ib][ie] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-}\r
-\r
-static void test_intdos(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig,\r
-  double* mat\r
-) {\r
-  int ib, ik, ne, ie;\r
-  double*** wght, e0[5], val0[2][5], val[2][5];\r
-\r
-  ne = 5;\r
-  wght = (double***)malloc(nk * sizeof(double**));\r
-  wght[0] = (double**)malloc(nk * nb * sizeof(double*));\r
-  wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;\r
-    }\r
-  }\r
-\r
-  for (ie = 0; ie < ne; ie++) {\r
-    e0[ie] = 0.2 * (double)(ie + 1);\r
-    val0[0][ie] = 4.0 * pi * pow(e0[ie], 5) / 5.0;\r
-    if (e0[ie] > 1.0 / sqrt(2.0))\r
-      val0[1][ie] = pi * pow(sqrt(-1.0 + 2.0 * pow(e0[ie], 2)), 5) / (5.0 * sqrt(2.0));\r
-    else\r
-      val0[1][ie] = 0.0;\r
-    e0[ie] = pow(e0[ie], 2) * 0.5;\r
-  }\r
-\r
-  intdos(ng, nk, nb, ne, bvec, eig, e0, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (ie = 0; ie < ne; ie++)\r
-      val[ib][ie] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (ie = 0; ie < ne; ie++)\r
-        val[ib][ie] += wght[ik][ib][ie] * mat[ik];\r
-\r
-  printf("# libtetrabz_intdos\n");\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (ie = 0; ie < ne; ie++)\r
-      printf("     %15.5e %15.5e\n", val0[ib][ie], val[ib][ie] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-}\r
-\r
-\r
-static void test_dblstep(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig1,\r
-  double** eig2,\r
-  double* mat\r
-) {\r
-  int ib, ik, jb;\r
-  double*** wght, val[2][2], ** eig12, ** eig22;\r
-\r
-  wght = (double***)malloc(nk * sizeof(double**));\r
-  eig12 = (double**)malloc(nk * sizeof(double*));\r
-  eig22 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double**)malloc(nk * nb * sizeof(double*));\r
-  eig12[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig22[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig12[ik] = eig12[0] + ik * nb;\r
-    eig22[ik] = eig22[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;\r
-    }\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++) {\r
-      eig12[ik][ib] = eig1[ik][ib] - 0.5;\r
-      eig22[ik][ib] = eig2[ik][ib] - 0.5;\r
-    }\r
-\r
-  dblstep(ng, nk, nb, bvec, eig12, eig22, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      val[ib][jb] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (jb = 0; jb < nb; jb++)\r
-        val[ib][jb] += wght[ik][ib][jb] * mat[ik];\r
-\r
-  printf("# libtetrabz_dblstep\n");\r
-  printf("     %15.5e %15.5e\n", 49.0 * pi / 320.0, val[0][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", 0.0, val[0][1] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi * (512.0 * sqrt(2.0) - 319.0) / 10240.0, val[1][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", 0.0, val[1][1] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig12[0]);\r
-  free(eig12);\r
-  free(eig22[0]);\r
-  free(eig22);\r
-}\r
-\r
-\r
-static void test_dbldelta(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig1,\r
-  double** eig2,\r
-  double* mat\r
-) {\r
-  int ib, ik, jb;\r
-  double*** wght, val[2][2], ** eig12, ** eig22;\r
-\r
-  wght = (double***)malloc(nk * sizeof(double**));\r
-  eig12 = (double**)malloc(nk * sizeof(double*));\r
-  eig22 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double**)malloc(nk * nb * sizeof(double*));\r
-  eig12[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig22[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig12[ik] = eig12[0] + ik * nb;\r
-    eig22[ik] = eig22[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;\r
-    }\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++) {\r
-      eig12[ik][ib] = eig1[ik][ib] - 0.5;\r
-      eig22[ik][ib] = eig2[ik][ib] - 0.5;\r
-    }\r
-\r
-  dbldelta(ng, nk, nb, bvec, eig12, eig22, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      val[ib][jb] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (jb = 0; jb < nb; jb++)\r
-        val[ib][jb] += wght[ik][ib][jb] * mat[ik];\r
-  printf("# libtetrabz_dbldelta\n");\r
-  printf("     %15.5e %15.5e\n", 2.0 * pi, val[0][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", 0.0, val[0][1] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi, val[1][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", 0.0, val[1][1] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig12[0]);\r
-  free(eig12);\r
-  free(eig22[0]);\r
-  free(eig22);\r
-}\r
-\r
-\r
-static void test_polstat(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig1,\r
-  double** eig2,\r
-  double* mat\r
-) {\r
-  int ib, ik, jb;\r
-  double*** wght, val[2][2], ** eig12, ** eig22;\r
-\r
-  wght = (double***)malloc(nk * sizeof(double**));\r
-  eig12 = (double**)malloc(nk * sizeof(double*));\r
-  eig22 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double**)malloc(nk * nb * sizeof(double*));\r
-  eig12[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig22[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig12[ik] = eig12[0] + ik * nb;\r
-    eig22[ik] = eig22[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;\r
-    }\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++) {\r
-      eig12[ik][ib] = eig1[ik][ib] - 0.5;\r
-      eig22[ik][ib] = eig2[ik][ib] - 0.5;\r
-    }\r
-\r
-  polstat(ng, nk, nb, bvec, eig12, eig22, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      val[ib][jb] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (jb = 0; jb < nb; jb++)\r
-        val[ib][jb] += wght[ik][ib][jb] * mat[ik];\r
-  printf("# libtetrabz_polstat\n");\r
-  printf("     %15.5e %15.5e\n", pi * (68.0 + 45.0 * log(3.0)) / 96.0, val[0][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi * 8.0 / 5.0, val[0][1] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi * (228.0 + 22.0 * sqrt(2.0) - 96.0 * log(2.0)\r
-    + 192.0 * log(4.0 + sqrt(2.0))\r
-    - 3.0 * log(1.0 + 2.0 * sqrt(2.0))) / 1536.0,\r
-    val[1][0] * vbz);\r
-  printf("     %15.5e %15.5e\n", pi * sqrt(8.0) / 5.0, val[1][1] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig12[0]);\r
-  free(eig12);\r
-  free(eig22[0]);\r
-  free(eig22);\r
-}\r
-\r
-\r
-static void test_fermigr(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig1,\r
-  double** eig2,\r
-  double* mat\r
-) {\r
-  int ib, ik, jb, ne, ie;\r
-  double**** wght, val[2][2][3], val0[2][2][3], ** eig12, ** eig22, e0[3];\r
-\r
-  ne = 3;\r
-\r
-  wght = (double****)malloc(nk * sizeof(double***));\r
-  eig12 = (double**)malloc(nk * sizeof(double*));\r
-  eig22 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double***)malloc(nk * nb * sizeof(double**));\r
-  eig12[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig22[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  wght[0][0] = (double**)malloc(nk * nb * nb * sizeof(double*));\r
-  wght[0][0][0] = (double*)malloc(nk * nb * nb * ne * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig12[ik] = eig12[0] + ik * nb;\r
-    eig22[ik] = eig22[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;\r
-      for (jb = 0; jb < nb; jb++) {\r
-        wght[ik][ib][jb] = wght[0][0][0] + ik * nb * nb * ne + ib * nb * ne + jb * ne;\r
-      }\r
-    }\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++) {\r
-      eig12[ik][ib] = eig1[ik][ib] - 0.5;\r
-      eig22[ik][ib] = eig2[ik][ib] - 0.5;\r
-    }\r
-\r
-  e0[0] = 1.0 / 3.0;\r
-  e0[1] = 2.0 / 3.0;\r
-  e0[2] = 1.0;\r
-  val0[0][0][0] = 4.0 * pi / 9.0;\r
-  val0[0][0][1] = 1295.0 * pi / 2592.0;\r
-  val0[0][0][2] = 15.0 * pi / 32.0;\r
-  val0[1][0][0] = 5183.0 * pi / 41472.0;\r
-  val0[1][0][1] = 4559.0 * pi / 41472.0;\r
-  val0[1][0][2] = 0.0;\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (ie = 0; ie < ne; ie++)\r
-      val0[ib][1][ie] = 0.0;\r
-\r
-  fermigr(ng, nk, nb, ne, bvec, eig12, eig22, e0, wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      for (ie = 0; ie < ne; ie++)\r
-        val[ib][jb][ie] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (jb = 0; jb < nb; jb++)\r
-        for (ie = 0; ie < ne; ie++)\r
-          val[ib][jb][ie] += wght[ik][ib][jb][ie] * mat[ik];\r
-\r
-  printf("# libtetrabz_fermigr\n");\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      for (ie = 0; ie < ne; ie++)\r
-        printf("     %15.5e %15.5e\n", val0[ib][jb][ie], val[ib][jb][ie] * vbz);\r
-  printf("\n");\r
-\r
-  free(wght[0][0][0]);\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig12[0]);\r
-  free(eig12);\r
-  free(eig22[0]);\r
-  free(eig22);\r
-}\r
-\r
-static void test_polcmplx(\r
-  int ng[3],\r
-  int nk,\r
-  int nb,\r
-  double bvec[3][3],\r
-  double vbz,\r
-  double** eig1,\r
-  double** eig2,\r
-  double* mat\r
-) {\r
-  int ib, ik, jb, ne, ie, ir;\r
-  double***** wght, val[2][2][3][2], val0[2][2][3][2], ** eig12, ** eig22, **e0;\r
-\r
-  ne = 3;\r
-\r
-  e0 = (double**)malloc(ne * sizeof(double*));\r
-  e0[0] = (double*)malloc(ne * 2 * sizeof(double));\r
-  for (ie = 0; ie < ne; ie++)\r
-    e0[ie] = e0[0] + ie * 2;\r
-\r
-  wght = (double*****)malloc(nk * sizeof(double****));\r
-  eig12 = (double**)malloc(nk * sizeof(double*));\r
-  eig22 = (double**)malloc(nk * sizeof(double*));\r
-  wght[0] = (double****)malloc(nk * nb * sizeof(double***));\r
-  eig12[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig22[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  wght[0][0] = (double***)malloc(nk * nb * nb * sizeof(double**));\r
-  wght[0][0][0] = (double**)malloc(nk * nb * nb * ne * sizeof(double*));\r
-  wght[0][0][0][0] = (double*)malloc(nk * nb * nb * ne * 2 * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    wght[ik] = wght[0] + ik * nb;\r
-    eig12[ik] = eig12[0] + ik * nb;\r
-    eig22[ik] = eig22[0] + ik * nb;\r
-    for (ib = 0; ib < nb; ib++) {\r
-      wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;\r
-      for (jb = 0; jb < nb; jb++) {\r
-        wght[ik][ib][jb] = wght[0][0][0] + ik * nb * nb * ne + ib * nb * ne + jb * ne;\r
-        for (ie = 0; ie < ne; ie++) {\r
-          wght[ik][ib][jb][ie] = wght[0][0][0][0] + ik * nb * nb * ne * 2\r
-            + ib * nb * ne * 2 + jb * ne * 2 + ie * 2;\r
-        }\r
-      }\r
-    }\r
-  }\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++) {\r
-      eig12[ik][ib] = eig1[ik][ib] - 0.5;\r
-      eig22[ik][ib] = eig2[ik][ib] - 0.5;\r
-    }\r
-\r
-  e0[0][0] = -2.0;\r
-  e0[0][1] = 1.0;\r
-  e0[1][0] = 0.0;\r
-  e0[1][1] = 2.0;\r
-  e0[2][0] = 1.0;\r
-  e0[2][1] = -0.5;\r
-  val0[0][0][0][0] = -0.838243341280338;\r
-  val0[0][0][0][1] = -0.734201894333234;\r
-  val0[0][0][1][0] = 0.270393588876530;\r
-  val0[0][0][1][1] = -0.771908416949610;\r
-  val0[0][0][2][0] = 0.970996830573510;\r
-  val0[0][0][2][1] = 0.302792326476720;\r
-  val0[1][0][0][0] = -0.130765724778920;\r
-  val0[1][0][0][1] = -0.087431218706638;\r
-  val0[1][0][1][0] = 0.030121954547245;\r
-  val0[1][0][1][1] = -0.135354254293510;\r
-  val0[1][0][2][0] = 0.178882244951203;\r
-  val0[1][0][2][1] = 0.064232167683425;\r
-  \r
-  for (ie = 0; ie < ne; ie++) {\r
-    for (ir = 0; ir < 2; ir++) {\r
-      val0[0][1][ie][0] = (8.0 * pi * (1.0 + 2.0 * e0[ie][0])) / (5.0 * (pow(1.0 + 2.0 * e0[ie][0], 2) + pow(2.0 * e0[ie][1], 2)));\r
-      val0[0][1][ie][1] = (- 8.0 * pi * 2.0 * e0[ie][1]) / (5.0 * (pow(1.0 + 2.0 * e0[ie][0], 2) + pow(2.0 * e0[ie][1], 2)));\r
-      val0[1][1][ie][0] = (sqrt(8.0) * pi * (1.0 + 4.0 * e0[ie][0])) / (5.0 * (pow(1.0 + 4.0 * e0[ie][0],2) + pow(4.0 * e0[ie][1], 2)));\r
-      val0[1][1][ie][1] = (- sqrt(8.0) * pi * 4.0 * e0[ie][1]) / (5.0 * (pow(1.0 + 4.0 * e0[ie][0], 2) + pow(4.0 * e0[ie][1], 2)));\r
-    }\r
-  }\r
-  \r
-  polcmplx(ng, nk, nb, ne, bvec, eig12, eig22, &e0[0], wght);\r
-\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      for (ie = 0; ie < ne; ie++)\r
-        for (ir = 0; ir < 2; ir++)\r
-          val[ib][jb][ie][ir] = 0.0;\r
-  for (ik = 0; ik < nk; ik++)\r
-    for (ib = 0; ib < nb; ib++)\r
-      for (jb = 0; jb < nb; jb++)\r
-        for (ie = 0; ie < ne; ie++)\r
-          for (ir = 0; ir < 2; ir++)\r
-            val[ib][jb][ie][ir] += wght[ik][ib][jb][ie][ir] * mat[ik];\r
-\r
-\r
-  printf("# libtetrabz_polcmplx\n");\r
-  for (ib = 0; ib < nb; ib++)\r
-    for (jb = 0; jb < nb; jb++)\r
-      for (ie = 0; ie < ne; ie++) {\r
-        printf("     %15.5e %15.5e\n", val0[ib][jb][ie][0], val[ib][jb][ie][0] * vbz);\r
-        printf("     %15.5e %15.5e\n", val0[ib][jb][ie][1], val[ib][jb][ie][1] * vbz);\r
-      }\r
-  printf("\n");\r
-  free(wght[0][0][0]);\r
-  free(wght[0][0]);\r
-  free(wght[0]);\r
-  free(wght);\r
-  free(eig12[0]);\r
-  free(eig12);\r
-  free(eig22[0]);\r
-  free(eig22); \r
-  free(e0[0]);\r
-  free(e0);\r
-}\r
-\r
-int main()\r
-{\r
-  int ng0, ng[3], nb, i3, j3, nk, ik, i0, i1, i2;\r
-  double bvec[3][3], vbz, ** eig1, ** eig2, * mat, kvec[3];\r
-\r
-  ng0 = 8;\r
-  for (i3 = 0; i3 < 3; i3++)    ng[i3] = ng0;\r
-  nk = ng[0] * ng[1] * ng[2];\r
-  nb = 2;\r
-  for (i3 = 0; i3 < 3; i3++)\r
-    for (j3 = 0; j3 < 3; j3++)\r
-      bvec[i3][j3] = 0.0;\r
-  for (i3 = 0; i3 < 3; i3++) bvec[i3][i3] = 3.0;\r
-  vbz = 27.0;/*abs(numpy.linalg.det(bvec));*/\r
-\r
-  mat = (double*)malloc(nk * sizeof(double));\r
-  eig1 = (double**)malloc(nk * sizeof(double*));\r
-  eig2 = (double**)malloc(nk * sizeof(double*));\r
-  eig1[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  eig2[0] = (double*)malloc(nk * nb * sizeof(double));\r
-  for (ik = 0; ik < nk; ik++) {\r
-    eig1[ik] = eig1[0] + ik * nb;\r
-    eig2[ik] = eig2[0] + ik * nb;\r
-  }\r
-\r
-  for (i0 = 0; i0 < ng[0]; i0++) {\r
-    for (i1 = 0; i1 < ng[1]; i1++) {\r
-      for (i2 = 0; i2 < ng[2]; i2++) {\r
-\r
-        kvec[0] = (double)i0 / (double)ng[0];\r
-        kvec[1] = (double)i1 / (double)ng[1];\r
-        kvec[2] = (double)i2 / (double)ng[2];\r
-        for (i3 = 0; i3 < 3; i3++)\r
-          if (kvec[i3] > 0.5) kvec[i3] += -1.0;\r
-        for (i3 = 0; i3 < 3; i3++) kvec[i3] = kvec[i3] * bvec[i3][i3];/*kvec.dot(bvec)*/\r
-\r
-        ik = i2 + i1 * ng[2] + i0 * ng[1] * ng[2];\r
-        eig1[ik][0] = 0.5 * (kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2]);\r
-        eig1[ik][1] = eig1[ik][0] + 0.25;\r
-\r
-        kvec[0] += 1.0;\r
-        eig2[ik][0] = 0.5 * (kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2]);\r
-        eig2[ik][1] = eig1[ik][0] + 0.5;\r
-      }\r
-    }\r
-  }\r
-\r
-  for (i0 = 0; i0 < ng[0]; i0++) {\r
-    for (i1 = 0; i1 < ng[1]; i1++) {\r
-      for (i2 = 0; i2 < ng[2]; i2++) {\r
-\r
-        kvec[0] = (double)i0 / (double)ng[0];\r
-        kvec[1] = (double)i1 / (double)ng[1];\r
-        kvec[2] = (double)i2 / (double)ng[2];\r
-        for (i3 = 0; i3 < 3; i3++)\r
-          if (kvec[i3] > 0.5) kvec[i3] += -1.0;\r
-        for (i3 = 0; i3 < 3; i3++) kvec[i3] = kvec[i3] * bvec[i3][i3];/*kvec.dot(bvec)*/\r
-\r
-        ik = i2 + i1 * ng[2] + i0 * ng[1] * ng[2];\r
-        mat[ik] = kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2];\r
-      }\r
-    }\r
-  }\r
-\r
-  printf("#          Ideal          Result\n");\r
-  \r
-  test_occ(ng, nk, nb, bvec, vbz, eig1, mat);\r
-  \r
-  test_fermieng(ng, nk, nb, bvec, vbz, eig1, mat);\r
-\r
-  test_dos(ng, nk, nb, bvec, vbz, eig1, mat);\r
-  \r
-  test_intdos(ng, nk, nb, bvec, vbz, eig1, mat);\r
-  \r
-  test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat);\r
-  \r
-  test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat);\r
-  \r
-  test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat);\r
-  \r
-  test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat);\r
-  \r
-  test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat);\r
-}\r
diff --git a/src/c/test.py b/src/c/test.py
deleted file mode 100644 (file)
index b368628..0000000
+++ /dev/null
@@ -1,479 +0,0 @@
-#\r
-# Copyright (c) 2014 Mitsuaki Kawamura\r
-#\r
-# Permission is hereby granted, free of charge, to any person obtaining a\r
-# copy of this software and associated documentation files (the\r
-# "Software"), to deal in the Software without restriction, including\r
-# without limitation the rights to use, copy, modify, merge, publish,\r
-# distribute, sublicense, and/or sell copies of the Software, and to\r
-# permit persons to whom the Software is furnished to do so, subject to\r
-# the following conditions:\r
-#\r
-# The above copyright notice and this permission notice shall be included\r
-# in all copies or substantial portions of the Software.\r
-#\r
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS\r
-# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
-# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.\r
-# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY\r
-# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,\r
-# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE\r
-# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
-#\r
-import math\r
-import libtetrabzc\r
-\r
-\r
-def test_occ(ng, nk, nb, bvec, vbz, eig, mat):\r
-    """\r
-\r
-    :return:\r
-    """\r
-    #\r
-    eig0 = eig.copy()\r
-    for ikb in range(nk * nb):\r
-        eig0[ikb] += - 0.5\r
-    #\r
-    wght = libtetrabzc.occ(ng[0], ng[1], ng[2], nk, nb,\r
-                          bvec[0][0], bvec[0][1], bvec[0][2],\r
-                          bvec[1][0], bvec[1][1], bvec[1][2],\r
-                          bvec[2][0], bvec[2][1], bvec[2][2], eig0)\r
-    #\r
-    val = [0.0, 0.0]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            val[ib] += wght[ik * nb + ib] * mat[ik]\r
-    print("# libtetrabz_occ")\r
-    print('     %15.5e %15.5e' % (4.0 * math.pi / 5.0, val[0] * vbz))\r
-    print('     %15.5e %15.5e' % (math.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))\r
-    print("")\r
-\r
-\r
-def test_fermieng(ng, nk, nb, bvec, vbz, eig, mat):\r
-    """\r
-\r
-    :return:\r
-    """\r
-    #\r
-    nelec = (4.0 * math.pi / 3.0 + math.sqrt(2.0) * math.pi / 3.0) / vbz\r
-\r
-    wght = libtetrabzc.fermieng(ng[0], ng[1], ng[2], nk, nb,\r
-                               bvec[0][0], bvec[0][1], bvec[0][2],\r
-                               bvec[1][0], bvec[1][1], bvec[1][2],\r
-                               bvec[2][0], bvec[2][1], bvec[2][2], eig, nelec)\r
-    ef = wght[nk * nb]\r
-    # iterations = wght[nk*nb + 1]\r
-    #\r
-    val = [0.0, 0.0]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            val[ib] += wght[ik * nb + ib] * mat[ik]\r
-    #\r
-    print("# libtetrabz_fermieng")\r
-    print("     %15.5e %15.5e" % (0.5, ef))\r
-    print("     %15.5e %15.5e" % (4.0 * math.pi / 5.0, val[0] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))\r
-    print("")\r
-\r
-\r
-def test_dos(ng, nk, nb, bvec, vbz, eig, mat):\r
-    """\r
-\r
-    :return:\r
-    """\r
-    ne = 5\r
-    #\r
-    e0 = [0.0, 0.0, 0.0, 0.0, 0.0]\r
-    val0 = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]\r
-    for ie in range(ne):\r
-        e0[ie] = 0.2 * (ie + 1)\r
-        val0[0][ie] = 4.0 * math.pi * e0[ie] ** 3\r
-        if e0[ie] > 1.0 / math.sqrt(2.0):\r
-            val0[1][ie] = math.sqrt(2.0) * math.pi * math.sqrt(-1.0 + 2.0 * e0[ie] ** 2) ** 3\r
-        else:\r
-            val0[1][ie] = 0.0\r
-    for ie in range(ne):\r
-        e0[ie] = e0[ie] ** 2 * 0.5\r
-    #\r
-    wght = libtetrabzc.dos(ng[0], ng[1], ng[2], nk, nb, ne,\r
-                          bvec[0][0], bvec[0][1], bvec[0][2],\r
-                          bvec[1][0], bvec[1][1], bvec[1][2],\r
-                          bvec[2][0], bvec[2][1], bvec[2][2], eig, e0)\r
-    #\r
-    val = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for ie in range(ne):\r
-                val[ib][ie] += wght[ik * nb * ne + ib * ne + ie] * mat[ik]\r
-    #\r
-    print("# libtetrabz_dos")\r
-    for ib in range(nb):\r
-        for ie in range(ne):\r
-            print("     %15.5e %15.5e" % (val0[ib][ie], val[ib][ie] * vbz))\r
-    print("")\r
-\r
-\r
-def test_intdos(ng, nk, nb, bvec, vbz, eig, mat):\r
-    """\r
-\r
-    :return:\r
-    """\r
-    ne = 5\r
-    e0 = [0.0, 0.0, 0.0, 0.0, 0.0]\r
-    val0 = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]\r
-    #\r
-    for ie in range(ne):\r
-        e0[ie] = 0.2 * (ie + 1)\r
-        val0[0][ie] = 4.0 * math.pi * e0[ie] ** 5 / 5.0\r
-        if e0[ie] > 1.0 / math.sqrt(2.0):\r
-            val0[1][ie] = math.pi * math.sqrt(-1.0 + 2.0 * e0[ie] ** 2) ** 5 / (5.0 * math.sqrt(2.0))\r
-        else:\r
-            val0[1][ie] = 0.0\r
-    for ie in range(ne):\r
-        e0[ie] = e0[ie] ** 2 * 0.5\r
-    #\r
-    wght = libtetrabzc.intdos(ng[0], ng[1], ng[2], nk, nb, ne,\r
-                             bvec[0][0], bvec[0][1], bvec[0][2],\r
-                             bvec[1][0], bvec[1][1], bvec[1][2],\r
-                             bvec[2][0], bvec[2][1], bvec[2][2], eig, e0)\r
-    #\r
-    val = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for ie in range(ne):\r
-                val[ib][ie] += wght[ik * nb * ne + ib * ne + ie] * mat[ik]\r
-    #\r
-    print("# libtetrabz_intdos")\r
-    for ib in range(nb):\r
-        for ie in range(ne):\r
-            print("     %15.5e %15.5e" % (val0[ib][ie], val[ib][ie] * vbz))\r
-    print("")\r
-\r
-\r
-def test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat):\r
-    """\r
-\r
-    :param bvec:\r
-    :param vbz:\r
-    :param eig1:\r
-    :param eig2:\r
-    :param mat:\r
-    :return:\r
-    """\r
-    #\r
-    eig10 = eig1.copy()\r
-    eig20 = eig2.copy()\r
-    for ikb in range(nk * nb):\r
-        eig10[ikb] += - 0.5\r
-        eig20[ikb] += - 0.5\r
-\r
-    wght = libtetrabzc.dblstep(ng[0], ng[1], ng[2], nk, nb,\r
-                              bvec[0][0], bvec[0][1], bvec[0][2],\r
-                              bvec[1][0], bvec[1][1], bvec[1][2],\r
-                              bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)\r
-    #\r
-    val = [[0.0, 0.0], [0.0, 0.0]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for jb in range(nb):\r
-                val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]\r
-    #\r
-    print("# libtetrabz_dblstep")\r
-    print("     %15.5e %15.5e" % (49.0 * math.pi / 320.0, val[0][0] * vbz))\r
-    print("     %15.5e %15.5e" % (0.0, val[0][1] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi * (512.0 * math.sqrt(2.0) - 319.0) / 10240.0, val[1][0] * vbz))\r
-    print("     %15.5e %15.5e" % (0.0, val[1][1] * vbz))\r
-    print("")\r
-\r
-\r
-def test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat):\r
-    """\r
-\r
-    :param bvec:\r
-    :param vbz:\r
-    :param eig1:\r
-    :param eig2:\r
-    :param mat:\r
-    :return:\r
-    """\r
-    #\r
-    eig10 = eig1.copy()\r
-    eig20 = eig2.copy()\r
-    for ikb in range(nk * nb):\r
-        eig10[ikb] += - 0.5\r
-        eig20[ikb] += - 0.5\r
-\r
-    wght = libtetrabzc.dbldelta(ng[0], ng[1], ng[2], nk, nb,\r
-                               bvec[0][0], bvec[0][1], bvec[0][2],\r
-                               bvec[1][0], bvec[1][1], bvec[1][2],\r
-                               bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)\r
-    #\r
-    val = [[0.0, 0.0], [0.0, 0.0]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for jb in range(nb):\r
-                val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]\r
-    #\r
-    print("# libtetrabz_dbldelta")\r
-    print("     %15.5e %15.5e" % (2.0 * math.pi, val[0][0] * vbz))\r
-    print("     %15.5e %15.5e" % (0.0, val[0][1] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi, val[1][0] * vbz))\r
-    print("     %15.5e %15.5e" % (0.0, val[1][1] * vbz))\r
-    print("")\r
-\r
-\r
-def test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat):\r
-    """\r
-\r
-    :param bvec:\r
-    :param vbz:\r
-    :param eig1:\r
-    :param eig2:\r
-    :param mat:\r
-    :return:\r
-    """\r
-    eig10 = eig1.copy()\r
-    eig20 = eig2.copy()\r
-    for ikb in range(nk * nb):\r
-        eig10[ikb] += - 0.5\r
-        eig20[ikb] += - 0.5\r
-\r
-    wght = libtetrabzc.polstat(ng[0], ng[1], ng[2], nk, nb,\r
-                              bvec[0][0], bvec[0][1], bvec[0][2],\r
-                              bvec[1][0], bvec[1][1], bvec[1][2],\r
-                              bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)\r
-    #\r
-    val = [[0.0, 0.0], [0.0, 0.0]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for jb in range(nb):\r
-                val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]\r
-    #\r
-    print("# libtetrabz_polstat")\r
-    print("     %15.5e %15.5e" % (math.pi * (68.0 + 45.0 * math.log(3.0)) / 96.0, val[0][0] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi * 8.0 / 5.0, val[0][1] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi * (228.0 + 22.0 * math.sqrt(2.0) - 96.0 * math.log(2.0)\r
-                                             + 192.0 * math.log(4.0 + math.sqrt(2.0))\r
-                                             - 3.0 * math.log(1.0 + 2.0 * math.sqrt(2.0))) / 1536.0,\r
-                                  val[1][0] * vbz))\r
-    print("     %15.5e %15.5e" % (math.pi * math.sqrt(8.0) / 5.0, val[1][1] * vbz))\r
-    print("")\r
-\r
-\r
-def test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat):\r
-    """\r
-\r
-    :param bvec:\r
-    :param vbz:\r
-    :param eig1:\r
-    :param eig2:\r
-    :param mat:\r
-    :return:\r
-    """\r
-    ne = 3\r
-    #\r
-    eig10 = eig1.copy()\r
-    eig20 = eig2.copy()\r
-    for ikb in range(nk * nb):\r
-        eig10[ikb] += - 0.5\r
-        eig20[ikb] += - 0.5\r
-\r
-    e0 = [0.0, 0.0, 0.0]\r
-    val0 = [[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]]\r
-    e0[0] = 1.0 / 3.0\r
-    e0[1] = 2.0 / 3.0\r
-    e0[2] = 1.0\r
-    val0[0][0][0] = 4.0 * math.pi / 9.0\r
-    val0[0][0][1] = 1295.0 * math.pi / 2592.0\r
-    val0[0][0][2] = 15.0 * math.pi / 32.0\r
-    val0[1][0][0] = 5183.0 * math.pi / 41472.0\r
-    val0[1][0][1] = 4559.0 * math.pi / 41472.0\r
-    val0[1][0][2] = 0.0\r
-    for ib in range(nb):\r
-        for ie in range(ne):\r
-            val0[ib][1][ie] = 0.0\r
-\r
-    wght = libtetrabzc.fermigr(ng[0], ng[1], ng[2], nk, nb, ne,\r
-                              bvec[0][0], bvec[0][1], bvec[0][2],\r
-                              bvec[1][0], bvec[1][1], bvec[1][2],\r
-                              bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20, e0)\r
-    #\r
-    val = [[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for jb in range(nb):\r
-                for ie in range(ne):\r
-                    val[ib][jb][ie] += wght[ik * nb * nb * ne + ib * nb * ne + jb * ne + ie] * mat[ik]\r
-    #\r
-    print("# libtetrabz_fermigr")\r
-    for ib in range(nb):\r
-        for jb in range(nb):\r
-            for ie in range(ne):\r
-                print("     %15.5e %15.5e" % (val0[ib][jb][ie], val[ib][jb][ie] * vbz))\r
-    print("")\r
-\r
-\r
-def test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat):\r
-    """\r
-\r
-    :param bvec:\r
-    :param vbz:\r
-    :param eig1:\r
-    :param eig2:\r
-    :param mat:\r
-    :return:\r
-    """\r
-    ne = 3\r
-\r
-    eig10 = eig1.copy()\r
-    eig20 = eig2.copy()\r
-    for ikb in range(nk * nb):\r
-        eig10[ikb] += - 0.5\r
-        eig20[ikb] += - 0.5\r
-\r
-    e0 = [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]\r
-    e0[0][0] = -2.0\r
-    e0[0][1] = 1.0\r
-    e0[1][0] = 0.0\r
-    e0[1][1] = 2.0\r
-    e0[2][0] = 1.0\r
-    e0[2][1] = - 0.5\r
-    e01 = [e0[0][0], e0[0][1], e0[1][0], e0[1][1], e0[2][0], e0[2][1]]\r
-\r
-    val0 = [[[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]],\r
-            [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]]]\r
-    val0[0][0][0][0] = -0.838243341280338\r
-    val0[0][0][0][1] = - 0.734201894333234\r
-    val0[0][0][1][0] = 0.270393588876530\r
-    val0[0][0][1][1] = - 0.771908416949610\r
-    val0[0][0][2][0] = 0.970996830573510\r
-    val0[0][0][2][1] = 0.302792326476720\r
-    val0[1][0][0][0] = -0.130765724778920\r
-    val0[1][0][0][1] = - 0.087431218706638\r
-    val0[1][0][1][0] = 0.030121954547245\r
-    val0[1][0][1][1] = - 0.135354254293510\r
-    val0[1][0][2][0] = 0.178882244951203\r
-    val0[1][0][2][1] = 0.064232167683425\r
-    for ie in range(3):\r
-        val0[0][1][ie][0] = (8.0 * math.pi * (1.0 + 2.0 * e0[ie][0])) \\r
-                            / (5.0 * ((1.0 + 2.0 * e0[ie][0]) ** 2 + (2.0 * e0[ie][1]) ** 2))\r
-        val0[0][1][ie][1] = (- 8.0 * math.pi * 2.0 * e0[ie][1]) \\r
-                            / (5.0 * ((1.0 + 2.0 * e0[ie][0]) ** 2 + (2.0 * e0[ie][1]) ** 2))\r
-        val0[1][1][ie][0] = (math.sqrt(8.0) * math.pi * (1.0 + 4.0 * e0[ie][0])) \\r
-                            / (5.0 * ((1.0 + 4.0 * e0[ie][0]) ** 2 + (4.0 * e0[ie][1]) ** 2))\r
-        val0[1][1][ie][1] = (- math.sqrt(8.0) * math.pi * 4.0 * e0[ie][1]) \\r
-                            / (5.0 * ((1.0 + 4.0 * e0[ie][0]) ** 2 + (4.0 * e0[ie][1]) ** 2))\r
-\r
-    wght = libtetrabzc.polcmplx(ng[0], ng[1], ng[2], nk, nb, ne,\r
-                               bvec[0][0], bvec[0][1], bvec[0][2],\r
-                               bvec[1][0], bvec[1][1], bvec[1][2],\r
-                               bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20, e01)\r
-    #\r
-    val = [[[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]],\r
-           [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]]]\r
-    for ik in range(nk):\r
-        for ib in range(nb):\r
-            for jb in range(nb):\r
-                for ie in range(ne):\r
-                    for ir in range(2):\r
-                        val[ib][jb][ie][ir] += wght[ik * nb * nb * ne * 2\r
-                                                    + ib * nb * ne * 2 + jb * ne * 2 + ie * 2 + ir] * mat[ik]\r
-    #\r
-    print("# libtetrabz_polcmplx")\r
-    for ib in range(nb):\r
-        for jb in range(nb):\r
-            for ie in range(ne):\r
-                print("     %15.5e %15.5e" % (val0[ib][jb][ie][0], val[ib][jb][ie][0] * vbz))\r
-                print("     %15.5e %15.5e" % (val0[ib][jb][ie][1], val[ib][jb][ie][1] * vbz))\r
-    print("")\r
-\r
-\r
-def test():\r
-    ng0 = 8\r
-    ng = [ng0, ng0, ng0]\r
-    nk = ng[0] * ng[1] * ng[2]\r
-    nb = 2\r
-    bvec = [[3.0, 0.0, 0.0],\r
-            [0.0, 3.0, 0.0],\r
-            [0.0, 0.0, 3.0]]\r
-    vbz = bvec[0][0] * bvec[1][1] * bvec[2][2]  # abs(numpy.linalg.det(bvec))\r
-\r
-    #\r
-    eig1 = []\r
-    eig2 = []\r
-    mat = []\r
-    #\r
-    kvec = [0.0, 0.0, 0.0]\r
-    for i0 in range(ng[0]):\r
-        for i1 in range(ng[1]):\r
-            for i2 in range(ng[2]):\r
-                kvec[0] = i0 / ng[0]\r
-                kvec[1] = i1 / ng[1]\r
-                kvec[2] = i2 / ng[2]\r
-                #\r
-                for ii in range(3):\r
-                    if kvec[ii] > 0.5:\r
-                        kvec[ii] += -1\r
-                # kvec[0:3] = kvec.dot(bvec)\r
-                for ii in range(3):\r
-                    kvec[ii] *= bvec[ii][ii]\r
-                #\r
-                eig01 = 0.0\r
-                for ii in range(3):\r
-                    eig01 += kvec[ii] ** 2\r
-                eig01 *= 0.5\r
-                eig1.append(eig01)\r
-                eig1.append(eig01 + 0.25)\r
-                #\r
-                kvec[0] += 1.0\r
-                eig02 = 0.0\r
-                for ii in range(3):\r
-                    eig02 += kvec[ii] ** 2\r
-                eig02 *= 0.5\r
-                eig2.append(eig02)\r
-                eig2.append(eig01 + 0.5)\r
-                #\r
-    #\r
-    for i0 in range(ng[1]):\r
-        for i1 in range(ng[1]):\r
-            for i2 in range(ng[2]):\r
-                kvec[0] = i0 / ng[0]\r
-                kvec[1] = i1 / ng[1]\r
-                kvec[2] = i2 / ng[2]\r
-                #\r
-                for ii in range(3):\r
-                    if kvec[ii] > 0.5:\r
-                        kvec[ii] += -1\r
-                # kvec[0:3] = kvec.dot(bvec)\r
-                for ii in range(3):\r
-                    kvec[ii] *= bvec[ii][ii]\r
-                #\r
-                mat0 = 0.0\r
-                for ii in range(3):\r
-                    mat0 += kvec[ii] ** 2\r
-                #\r
-                mat.append(mat0)\r
-    #\r
-    print("#          Ideal          Result")\r
-    #\r
-    test_occ(ng, nk, nb, bvec, vbz, eig1, mat)\r
-    #\r
-    test_fermieng(ng, nk, nb, bvec, vbz, eig1, mat)\r
-    #\r
-    test_dos(ng, nk, nb, bvec, vbz, eig1, mat)\r
-    #\r
-    test_intdos(ng, nk, nb, bvec, vbz, eig1, mat)\r
-    #\r
-    test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat)\r
-    #\r
-    test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat)\r
-    #\r
-    test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat)\r
-    #\r
-    test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat)\r
-    #\r
-    test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat)\r
-    #\r
-\r
-\r
-test()\r