• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revision92ff29626bfd6f081ed3ead65e124d23083a36df (tree)
Zeit2022-03-18 01:58:11
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Ändern Zusammenfassung

Diff

--- a/src/c/libtetrabz_common.c
+++ b/src/c/libtetrabz_common.c
@@ -20,12 +20,12 @@
2020 // TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
2121 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 */
23-
23+#include <stdio.h>
2424
2525 void libtetrabz_initialize(
26- int *ng,
27- double **bvec,
28- double **wlsm,
26+ int ng[3],
27+ double bvec[3][3],
28+ double wlsm[20][4],
2929 int **ikv
3030 )
3131 {
@@ -157,6 +157,7 @@ void libtetrabz_initialize(
157157 ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
158158 ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
159159 ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
160+ for (i3 = 0; i3 < 3; i3++) if (ikv0[i3] < 0) ikv0[i3] += ng[i3];
160161 ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
161162 }
162163 nt += 1;
@@ -171,7 +172,7 @@ void libtetrabz_tsmall_a1(
171172 double *e,
172173 double e0,
173174 double *v,
174- double **tsmall
175+ double tsmall[4][4]
175176 ) {
176177 /*
177178 Cut small tetrahedron A1
@@ -206,7 +207,7 @@ void libtetrabz_tsmall_b1(
206207 double *e,
207208 double e0,
208209 double *v,
209- double **tsmall
210+ double tsmall[4][4]
210211 )
211212 {
212213 /*
@@ -244,7 +245,7 @@ void libtetrabz_tsmall_b2(
244245 double *e,
245246 double e0,
246247 double *v,
247- double **tsmall
248+ double tsmall[4][4]
248249 )
249250 {
250251 /*
@@ -280,7 +281,7 @@ void libtetrabz_tsmall_b3(
280281 double *e,
281282 double e0,
282283 double *v,
283- double **tsmall
284+ double tsmall[4][4]
284285 )
285286 {
286287 /*
@@ -318,7 +319,7 @@ void libtetrabz_tsmall_c1(
318319 double *e,
319320 double e0,
320321 double *v,
321- double **tsmall
322+ double tsmall[4][4]
322323 )
323324 {
324325 /*
@@ -354,7 +355,7 @@ void libtetrabz_tsmall_c2(
354355 double *e,
355356 double e0,
356357 double *v,
357- double **tsmall
358+ double tsmall[4][4]
358359 )
359360 {
360361 /*
@@ -390,7 +391,7 @@ void libtetrabz_tsmall_c3(
390391 double *e,
391392 double e0,
392393 double *v,
393- double **tsmall
394+ double tsmall[4][4]
394395 )
395396 {
396397 /*
@@ -428,7 +429,7 @@ void libtetrabz_triangle_a1(
428429 double *e,
429430 double e0,
430431 double *v,
431- double **tsmall
432+ double tsmall[4][3]
432433 )
433434 {
434435 /*
@@ -461,7 +462,7 @@ void libtetrabz_triangle_b1(
461462 double *e,
462463 double e0,
463464 double *v,
464- double **tsmall
465+ double tsmall[4][3]
465466 ) {
466467 /*
467468 Cut triangle B1
@@ -494,7 +495,7 @@ void libtetrabz_triangle_b2(
494495 double *e,
495496 double e0,
496497 double *v,
497- double **tsmall
498+ double tsmall[4][3]
498499 )
499500 {
500501 /*
@@ -529,7 +530,7 @@ void libtetrabz_triangle_c1(
529530 double *e,
530531 double e0,
531532 double *v,
532- double **tsmall
533+ double tsmall[4][3]
533534 )
534535 {
535536 /*
@@ -562,20 +563,30 @@ void eig_sort(
562563 int *swap //!< [out] Order of index (sorted)
563564 )
564565 {
565- int i, j, k;
566+ int i, j, k, min_loc;
567+ double min_val;
566568
567569 for (i = 0; i < n; ++i) swap[i] = i;
568570
569571 for (i = 0; i < n - 1; ++i) {
572+ min_val = key[i];
573+ min_loc = i;
570574 for (j = i + 1; j < n; ++j) {
571- if (key[swap[j]] < key[swap[i]]) {
572- /*
573- Swap
574- */
575- k = swap[j];
576- swap[j] = swap[i];
577- swap[i] = k;
575+ if (min_val > key[j]) {
576+ min_val = key[j];
577+ min_loc = j;
578578 }
579- }/*for (j = i + 1; j < n; ++j)*/
579+ }
580+ if (key[i] > min_val) {
581+ /*
582+ Swap
583+ */
584+ key[min_loc] = key[i];
585+ key[i] = min_val;
586+
587+ k = swap[min_loc];
588+ swap[min_loc] = swap[i];
589+ swap[i] = k;
590+ }
580591 }/*for (i = 0; i < n - 1; ++i)*/
581592 }/*eig_sort*/
--- a/src/c/libtetrabz_common.h
+++ b/src/c/libtetrabz_common.h
@@ -21,7 +21,7 @@
2121 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 */
2323
24-void libtetrabz_initialize(int *ng, double **bvec, double wlsm[4][20], int **ikv);
24+void libtetrabz_initialize(int ng[3], double bvec[3][3], double wlsm[20][4], int** ikv);
2525 void libtetrabz_tsmall_a1(double *e, double e0, double *v, double tsmall[4][4]);
2626 void libtetrabz_tsmall_b1(double *e, double e0, double *v, double tsmall[4][4]);
2727 void libtetrabz_tsmall_b2(double *e, double e0, double *v, double tsmall[4][4]);
@@ -29,8 +29,8 @@ void libtetrabz_tsmall_b3(double *e, double e0, double *v, double tsmall[4][4]);
2929 void libtetrabz_tsmall_c1(double *e, double e0, double *v, double tsmall[4][4]);
3030 void libtetrabz_tsmall_c2(double *e, double e0, double *v, double tsmall[4][4]);
3131 void libtetrabz_tsmall_c3(double *e, double e0, double *v, double tsmall[4][4]);
32-void libtetrabz_triangle_a1(double *e, double e0, double *v, double tsmall[3][4]);
33-void libtetrabz_triangle_b1(double *e, double e0, double *v, double tsmall[3][4]);
34-void libtetrabz_triangle_b2(double *e, double e0, double *v, double tsmall[3][4]);
35-void libtetrabz_triangle_c1(double *e, double e0, double *v, double tsmall[3][4]);
32+void libtetrabz_triangle_a1(double *e, double e0, double *v, double tsmall[4][3]);
33+void libtetrabz_triangle_b1(double *e, double e0, double *v, double tsmall[4][3]);
34+void libtetrabz_triangle_b2(double *e, double e0, double *v, double tsmall[4][3]);
35+void libtetrabz_triangle_c1(double *e, double e0, double *v, double tsmall[4][3]);
3636 void eig_sort(int n, double *key, int *swap);
--- a/src/c/libtetrabz_dbldelta.c
+++ b/src/c/libtetrabz_dbldelta.c
@@ -74,10 +74,10 @@ void libtetrabz_dbldelta2(
7474 }
7575
7676 void dbldelta(
77- int* ng,
77+ int ng[3],
7878 int nk,
7979 int nb,
80- double** bvec,
80+ double bvec[3][3],
8181 double** eig1,
8282 double** eig2,
8383 double*** wght
@@ -86,7 +86,7 @@ void dbldelta(
8686 /*
8787 Main SUBROUTINE for Delta(E1) * Delta(E2)
8888 */
89- int it, ik, ib, ii, jj, jb, ** ikv, indx[4];
89+ int it, ik, ib, i20, i4, j3, jb, ** ikv, indx[4];
9090 double wlsm[20][4], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[4][3], thr;
9191
9292 thr = 1.0e-10;
@@ -101,9 +101,9 @@ void dbldelta(
101101 ej1 = (double**)malloc(4 * sizeof(double*));
102102 ei1[0] = (double*)malloc(4 * nb * sizeof(double));
103103 ej1[0] = (double*)malloc(4 * nb * sizeof(double));
104- for (ii = 0; ii < 4; ii++) {
105- ei1[ii] = ei1[0] + ii * nb;
106- ej1[ii] = ej1[0] + ii * nb;
104+ for (i4 = 0; i4 < 4; i4++) {
105+ ei1[i4] = ei1[0] + i4 * nb;
106+ ej1[i4] = ej1[0] + i4 * nb;
107107 }
108108
109109 w1 = (double***)malloc(nb * sizeof(double**));
@@ -111,8 +111,8 @@ void dbldelta(
111111 w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
112112 for (ib = 0; ib < nb; ib++) {
113113 w1[ib] = w1[0] + ib * 4;
114- for (ii = 0; ii < 4; ii++) {
115- w1[ib][ii] = w1[0][0] + ib * 4 * nb + ii * nb;
114+ for (i4 = 0; i4 < 4; i4++) {
115+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
116116 }
117117 }
118118
@@ -137,27 +137,27 @@ void dbldelta(
137137
138138 for (it = 0; it < 6 * nk; it++) {
139139
140- for (ii = 0; ii < 4; ii++)
140+ for (i4 = 0; i4 < 4; i4++)
141141 for (ib = 0; ib < nb; ib++) {
142- ei1[ii][ib] = 0.0;
143- ej1[ii][ib] = 0.0;
142+ ei1[i4][ib] = 0.0;
143+ ej1[i4][ib] = 0.0;
144144 }
145- for (jj = 0; jj < 20; jj++) {
146- for (ii = 0; ii < 4; ii++) {
145+ for (i20 = 0; i20 < 20; i20++) {
146+ for (i4 = 0; i4 < 4; i4++) {
147147 for (ib = 0; ib < nb; ib++) {
148- ei1[ii][ib] += eig1[ikv[it][jj]][ib] * wlsm[jj][ii];
149- ej1[ii][ib] += eig2[ikv[it][jj]][ib] * wlsm[jj][ii];
148+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
149+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
150150 }
151151 }
152152 }
153153
154154 for (ib = 0; ib < nb; ib++) {
155155
156- for (ii = 0; ii < 4; ii++)
156+ for (i4 = 0; i4 < 4; i4++)
157157 for (jb = 0; jb < nb; jb++)
158- w1[ib][ii][jb] = 0.0;
158+ w1[ib][i4][jb] = 0.0;
159159
160- for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
160+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
161161 eig_sort(4, e, indx);
162162
163163 if (e[0] < 0.0 && e[0] <= e[1]) {
@@ -165,17 +165,17 @@ void dbldelta(
165165 libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
166166
167167 if (v > thr) {
168- for (jj = 0; jj < 3; jj++)
169- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
170- for (ii = 0; ii < 4; ii++)
171- for (jj = 0; jj < 3; jj++)
168+ for (j3 = 0; j3 < 3; j3++)
169+ for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
170+ for (i4 = 0; i4 < 4; i4++)
171+ for (j3 = 0; j3 < 3; j3++)
172172 for (jb = 0; jb < nb; jb++)
173- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
173+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
174174 libtetrabz_dbldelta2(nb, ej2, w2);
175- for (ii = 0; ii < 4; ii++)
175+ for (i4 = 0; i4 < 4; i4++)
176176 for (jb = 0; jb < nb; jb++)
177- for (jj = 0; jj < 3; jj++)
178- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
177+ for (j3 = 0; j3 < 3; j3++)
178+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
179179 }
180180 }
181181 else if (e[1] < 0.0 && 0.0 <= e[2]) {
@@ -183,33 +183,33 @@ void dbldelta(
183183 libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
184184
185185 if (v > thr) {
186- for (jj = 0; jj < 3; jj++)
187- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
188- for (ii = 0; ii < 4; ii++)
189- for (jj = 0; jj < 3; jj++)
186+ for (j3 = 0; j3 < 3; j3++)
187+ for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
188+ for (i4 = 0; i4 < 4; i4++)
189+ for (j3 = 0; j3 < 3; j3++)
190190 for (jb = 0; jb < nb; jb++)
191- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
191+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
192192 libtetrabz_dbldelta2(nb, ej2, w2);
193- for (ii = 0; ii < 4; ii++)
193+ for (i4 = 0; i4 < 4; i4++)
194194 for (jb = 0; jb < nb; jb++)
195- for (jj = 0; jj < 3; jj++)
196- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
195+ for (j3 = 0; j3 < 3; j3++)
196+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
197197 }
198198
199199 libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
200200
201201 if (v > thr) {
202- for (jj = 0; jj < 3; jj++)
203- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
204- for (ii = 0; ii < 4; ii++)
205- for (jj = 0; jj < 3; jj++)
202+ for (j3 = 0; j3 < 3; j3++)
203+ for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
204+ for (i4 = 0; i4 < 4; i4++)
205+ for (j3 = 0; j3 < 3; j3++)
206206 for (jb = 0; jb < nb; jb++)
207- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
207+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
208208 libtetrabz_dbldelta2(nb, ej2, w2);
209- for (ii = 0; ii < 4; ii++)
209+ for (i4 = 0; i4 < 4; i4++)
210210 for (jb = 0; jb < nb; jb++)
211- for (jj = 0; jj < 3; jj++)
212- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
211+ for (j3 = 0; j3 < 3; j3++)
212+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
213213 }
214214 }
215215 else if (e[2] < 0.0 && 0.0 < e[3]) {
@@ -217,28 +217,28 @@ void dbldelta(
217217 libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
218218
219219 if (v > thr) {
220- for (jj = 0; jj < 3; jj++)
221- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
222- for (ii = 0; ii < 4; ii++)
223- for (jj = 0; jj < 3; jj++)
220+ for (j3 = 0; j3 < 3; j3++)
221+ for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
222+ for (i4 = 0; i4 < 4; i4++)
223+ for (j3 = 0; j3 < 3; j3++)
224224 for (jb = 0; jb < nb; jb++)
225- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
225+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
226226 libtetrabz_dbldelta2(nb, ej2, w2);
227- for (ii = 0; ii < 4; ii++)
227+ for (i4 = 0; i4 < 4; i4++)
228228 for (jb = 0; jb < nb; jb++)
229- for (jj = 0; jj < 3; jj++)
230- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
229+ for (j3 = 0; j3 < 3; j3++)
230+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
231231 }
232232 }
233233 else {
234234 continue;
235235 }
236236 }
237- for (ii = 0; ii < 20; ii++)
237+ for (i20 = 0; i20 < 20; i20++)
238238 for (ib = 0; ib < nb; ib++)
239- for (jj = 0; jj < 4; jj++)
239+ for (i4 = 0; i4 < 4; i4++)
240240 for (jb = 0; jb < nb; jb++)
241- wght[ikv[it][ii]][ib][jb] += wlsm[ii][jj] * w1[ib][jj][jb];
241+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
242242 }
243243 for (ik = 0; ik < nk; ik++)
244244 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_dblstep.c
+++ b/src/c/libtetrabz_dblstep.c
@@ -26,7 +26,7 @@
2626
2727 void libtetrabz_dblstep2(
2828 int nb,
29- double *ei,
29+ double ei[4],
3030 double **ej,
3131 double **w1
3232 ) {
@@ -40,10 +40,10 @@ void libtetrabz_dblstep2(
4040
4141 for (ib = 0; ib < nb; ib++) {
4242
43- for (ii = 0; ii < 3; ii++)
43+ for (ii = 0; ii < 4; ii++)
4444 w1[ib][ii] = 0.0;
4545
46- for (ii = 0; ii < 3; ii++) e[ii] = - ei[ii] + ej[ib][ii];
46+ for (ii = 0; ii < 4; ii++) e[ii] = - ei[ii] + ej[ib][ii];
4747 eig_sort(4, e, indx);
4848
4949 if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
@@ -52,60 +52,60 @@ void libtetrabz_dblstep2(
5252 */
5353 v = 0.5;
5454 for (ii = 0; ii < 4; ii++)
55- w1[ib][ii] = v * 0.25;
55+ w1[ib][ii] += v * 0.25;
5656 }
5757 else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
5858 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
5959 for (ii = 0; ii < 4; ii++)
6060 for (jj = 0; jj < 4; jj++)
61- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
61+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
6262 }
6363 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
6464
6565 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
6666 for (ii = 0; ii < 4; ii++)
6767 for (jj = 0; jj < 4; jj++)
68- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
68+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
6969
7070 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
7171 for (ii = 0; ii < 4; ii++)
7272 for (jj = 0; jj < 4; jj++)
73- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
73+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
7474
7575 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
7676 for (ii = 0; ii < 4; ii++)
7777 for (jj = 0; jj < 4; jj++)
78- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
78+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
7979 }
8080 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
8181
8282 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
8383 for (ii = 0; ii < 4; ii++)
8484 for (jj = 0; jj < 4; jj++)
85- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
85+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
8686
8787 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
8888 for (ii = 0; ii < 4; ii++)
8989 for (jj = 0; jj < 4; jj++)
90- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
90+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
9191
9292 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
9393 for (ii = 0; ii < 4; ii++)
9494 for (jj = 0; jj < 4; jj++)
95- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
95+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
9696 }
9797 else if (e[3] <= 0.0) {
9898 for (ii = 0; ii < 4; ii++)
99- w1[ib][ii] = 0.25;
99+ w1[ib][ii] += 0.25;
100100 }
101101 }
102102 }
103103
104104 void dblstep(
105- int *ng,
105+ int ng[3],
106106 int nk,
107107 int nb,
108- double **bvec,
108+ double bvec[3][3],
109109 double **eig1,
110110 double **eig2,
111111 double ***wght
@@ -199,7 +199,7 @@ void dblstep(
199199 }
200200 for (ii = 0; ii < 4; ii++)
201201 for (jj = 0; jj < 4; jj++) {
202- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
202+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
203203 for (jb = 0; jb < nb; jb++)
204204 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
205205 }
@@ -221,7 +221,7 @@ void dblstep(
221221 }
222222 for (ii = 0; ii < 4; ii++)
223223 for (jj = 0; jj < 4; jj++) {
224- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
224+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
225225 for (jb = 0; jb < nb; jb++)
226226 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
227227 }
@@ -241,7 +241,7 @@ void dblstep(
241241 }
242242 for (ii = 0; ii < 4; ii++)
243243 for (jj = 0; jj < 4; jj++) {
244- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
244+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
245245 for (jb = 0; jb < nb; jb++)
246246 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
247247 }
@@ -261,7 +261,7 @@ void dblstep(
261261 }
262262 for (ii = 0; ii < 4; ii++)
263263 for (jj = 0; jj < 4; jj++) {
264- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
264+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
265265 for (jb = 0; jb < nb; jb++)
266266 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
267267 }
@@ -283,7 +283,7 @@ void dblstep(
283283 }
284284 for (ii = 0; ii < 4; ii++)
285285 for (jj = 0; jj < 4; jj++) {
286- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
286+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
287287 for (jb = 0; jb < nb; jb++)
288288 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
289289 }
@@ -303,7 +303,7 @@ void dblstep(
303303 }
304304 for (ii = 0; ii < 4; ii++)
305305 for (jj = 0; jj < 4; jj++) {
306- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
306+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
307307 for (jb = 0; jb < nb; jb++)
308308 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
309309 }
@@ -323,7 +323,7 @@ void dblstep(
323323 }
324324 for (ii = 0; ii < 4; ii++)
325325 for (jj = 0; jj < 4; jj++) {
326- ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
326+ ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
327327 for (jb = 0; jb < nb; jb++)
328328 ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
329329 }
@@ -336,7 +336,7 @@ void dblstep(
336336 }
337337 else if (e[3] <= 0.0) {
338338 for (ii = 0; ii < 4; ii++) {
339- ei2[ii] = ej1[ii][ib];
339+ ei2[ii] = ei1[ii][ib];
340340 for (jb = 0; jb < nb; jb++)
341341 ej2[jb][ii] = ej1[ii][jb];
342342 }
@@ -346,14 +346,14 @@ void dblstep(
346346 w1[ib][ii][jb] += w2[jb][ii];
347347 }
348348 else {
349- continue;
349+ continue;
350350 }
351351 }
352352 for (ii = 0; ii < 20; ii++)
353353 for (jj = 0; jj < 4; jj++)
354354 for (ib = 0; ib < nb; ib++)
355355 for (jb = 0; jb < nb; jb++)
356- wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
356+ wght[ikv[it][ii]][ib][jb] += w1[ib][jj][jb] * wlsm[ii][jj];
357357 }
358358 for (ik = 0; ik < nk; ik++)
359359 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_dos.c
+++ b/src/c/libtetrabz_dos.c
@@ -24,11 +24,11 @@
2424 #include <stdlib.h>
2525
2626 void dos(
27- int* ng,
27+ int ng[3],
2828 int nk,
2929 int nb,
3030 int ne,
31- double** bvec,
31+ double bvec[3][3],
3232 double** eig,
3333 double* e0,
3434 double*** wght) {
@@ -54,7 +54,7 @@ void dos(
5454 w1[0] = (double**)malloc(nb * ne * sizeof(double*));
5555 w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
5656 for (ib = 0; ib < nb; ib++) {
57- w1[ii] = w1[0] + ib * ne;
57+ w1[ib] = w1[0] + ib * ne;
5858 for (ie = 0; ie < ne; ie++) {
5959 w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
6060 }
@@ -77,13 +77,12 @@ void dos(
7777 for (ib = 0; ib < nb; ib++)
7878 ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
7979
80- for (ib = 0; ib < nb; ib++)
80+ for (ib = 0; ib < nb; ib++) {
81+
8182 for (ie = 0; ie < ne; ie++)
8283 for (ii = 0; ii < 4; ii++)
8384 w1[ib][ie][ii] = 0.0;
8485
85- for (ib = 0; ib < nb; ib++) {
86-
8786 for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
8887 eig_sort(4, e, indx);
8988
@@ -94,7 +93,7 @@ void dos(
9493 libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
9594 for (ii = 0; ii < 4; ii++)
9695 for (jj = 0; jj < 3; jj++)
97- w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
96+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
9897
9998 }
10099 else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
@@ -102,30 +101,30 @@ void dos(
102101 libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
103102 for (ii = 0; ii < 4; ii++)
104103 for (jj = 0; jj < 3; jj++)
105- w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
104+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
106105
107106 libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
108107 for (ii = 0; ii < 4; ii++)
109108 for (jj = 0; jj < 3; jj++)
110- w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
109+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
111110 }
112111 else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
113112
114113 libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
115114 for (ii = 0; ii < 4; ii++)
116115 for (jj = 0; jj < 3; jj++)
117- w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
116+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
118117 }
119118 else {
120119 continue;
121120 }
122- for (ii = 0; ii < 20; ii++)
123- for (ib = 0; ib < nb; ib++)
124- for (ie = 0; ie < ne; ie++)
125- for (jj = 0; jj < 4; jj++)
126- wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
127121 }
128122 }
123+ for (ii = 0; ii < 20; ii++)
124+ for (ib = 0; ib < nb; ib++)
125+ for (ie = 0; ie < ne; ie++)
126+ for (jj = 0; jj < 4; jj++)
127+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
129128 }
130129 for (ik = 0; ik < nk; ik++)
131130 for (ib = 0; ib < nb; ib++)
@@ -142,11 +141,11 @@ void dos(
142141 }
143142
144143 void intdos(
145- int* ng,
144+ int ng[3],
146145 int nk,
147146 int nb,
148147 int ne,
149- double** bvec,
148+ double bvec[3][3],
150149 double** eig,
151150 double* e0,
152151 double*** wght
@@ -174,7 +173,7 @@ void intdos(
174173 w1[0] = (double**)malloc(nb * ne * sizeof(double*));
175174 w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
176175 for (ib = 0; ib < nb; ib++) {
177- w1[ii] = w1[0] + ib * ne;
176+ w1[ib] = w1[0] + ib * ne;
178177 for (ie = 0; ie < ne; ie++) {
179178 w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
180179 }
@@ -213,24 +212,24 @@ void intdos(
213212 libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
214213 for (ii = 0; ii < 4; ii++)
215214 for (jj = 0; jj < 4; jj++)
216- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
215+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
217216 }
218217 else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
219218
220219 libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
221220 for (ii = 0; ii < 4; ii++)
222221 for (jj = 0; jj < 4; jj++)
223- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
222+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
224223
225224 libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
226225 for (ii = 0; ii < 4; ii++)
227226 for (jj = 0; jj < 4; jj++)
228- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
227+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
229228
230229 libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
231230 for (ii = 0; ii < 4; ii++)
232231 for (jj = 0; jj < 4; jj++)
233- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
232+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
234233
235234 }
236235 else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
@@ -238,17 +237,17 @@ void intdos(
238237 libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
239238 for (ii = 0; ii < 4; ii++)
240239 for (jj = 0; jj < 4; jj++)
241- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
240+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
242241
243242 libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
244243 for (ii = 0; ii < 4; ii++)
245244 for (jj = 0; jj < 4; jj++)
246- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
245+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
247246
248247 libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
249248 for (ii = 0; ii < 4; ii++)
250249 for (jj = 0; jj < 4; jj++)
251- w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
250+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
252251
253252 }
254253 else if (e[3] <= e0[ie]) {
--- a/src/c/libtetrabz_fermigr.c
+++ b/src/c/libtetrabz_fermigr.c
@@ -93,7 +93,7 @@ void libtetrabz_fermigr2(
9393 for (ie = 0; ie < ne; ie++)
9494 w1[ib][i4][ie] = 0.0;
9595
96- for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[i4][ib];
96+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
9797 eig_sort(4, e, indx);
9898
9999 if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
@@ -202,7 +202,7 @@ void libtetrabz_fermigr2(
202202 }
203203 else if (e[3] <= 0.0) {
204204 for (i4 = 0; i4 < 4; i4++) {
205- de[i4] = tsmall[i4][j4] * (ej1[i4][ib] - ei1[i4]);
205+ de[i4] = ej1[ib][i4] - ei1[i4];
206206 }
207207 libtetrabz_fermigr3(ne, e0, de, w2);
208208 for (i4 = 0; i4 < 4; i4++)
@@ -216,11 +216,11 @@ void libtetrabz_fermigr2(
216216 }
217217
218218 void fermigr(
219- int *ng,
219+ int ng[3],
220220 int nk,
221221 int nb,
222222 int ne,
223- double **bvec,
223+ double bvec[3][3],
224224 double **eig1,
225225 double **eig2,
226226 double *e0,
--- a/src/c/libtetrabz_occ.c
+++ b/src/c/libtetrabz_occ.c
@@ -22,13 +22,14 @@
2222 //
2323 #include <math.h>
2424 #include <stdlib.h>
25+#include <stdio.h>
2526 #include "libtetrabz_common.h"
2627
2728 void occ(
28- int* ng,
29+ int ng[3],
2930 int nk,
3031 int nb,
31- double** bvec,
32+ double bvec[3][3],
3233 double** eig,
3334 double** wght
3435 ) {
@@ -72,12 +73,11 @@ void occ(
7273 for (ib = 0; ib < nb; ib++)
7374 ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
7475
75- for (ib = 0; ib < nb; ib++)
76+ for (ib = 0; ib < nb; ib++) {
77+
7678 for (ii = 0; ii < 4; ii++)
7779 w1[ib][ii] = 0.0;
7880
79- for (ib = 0; ib < nb; ib++) {
80-
8181 for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
8282 eig_sort(4, e, indx);
8383
@@ -85,41 +85,41 @@ void occ(
8585 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
8686 for (ii = 0; ii < 4; ii++)
8787 for (jj = 0; jj < 4; jj++)
88- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
88+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
8989 }
9090 else if (e[1] <= 0.0 && 0.0 < e[2]) {
9191
9292 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
9393 for (ii = 0; ii < 4; ii++)
9494 for (jj = 0; jj < 4; jj++)
95- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
95+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
9696
9797 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
9898 for (ii = 0; ii < 4; ii++)
9999 for (jj = 0; jj < 4; jj++)
100- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
100+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
101101
102102 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
103103 for (ii = 0; ii < 4; ii++)
104104 for (jj = 0; jj < 4; jj++)
105- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
105+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
106106 }
107107 else if (e[2] <= 0.0 && 0.0 < e[3]) {
108108
109109 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
110110 for (ii = 0; ii < 4; ii++)
111111 for (jj = 0; jj < 4; jj++)
112- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
112+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
113113
114114 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
115115 for (ii = 0; ii < 4; ii++)
116116 for (jj = 0; jj < 4; jj++)
117- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
117+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
118118
119119 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
120120 for (ii = 0; ii < 4; ii++)
121121 for (jj = 0; jj < 4; jj++)
122- w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
122+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
123123 }
124124 else if (e[3] <= 0.0) {
125125 for (ii = 0; ii < 4; ii++)
@@ -131,8 +131,9 @@ void occ(
131131 }
132132 for (ii = 0; ii < 20; ii++)
133133 for (ib = 0; ib < nb; ib++)
134- for (jj = 0; jj < 4; jj++)
134+ for (jj = 0; jj < 4; jj++) {
135135 wght[ikv[it][ii]][ib] += wlsm[ii][jj] * w1[ib][jj];
136+ }
136137 }
137138 for (ik = 0; ik < nk; ik++)
138139 for (ib = 0; ib < nb; ib++)
@@ -147,10 +148,10 @@ void occ(
147148 }
148149
149150 void fermieng(
150- int *ng,
151+ int ng[3],
151152 int nk,
152153 int nb,
153- double **bvec,
154+ double bvec[3][3],
154155 double **eig,
155156 double **wght,
156157 double nelec,
--- a/src/c/libtetrabz_polcmplx.c
+++ b/src/c/libtetrabz_polcmplx.c
@@ -23,6 +23,7 @@
2323 #include "libtetrabz_common.h"
2424 #include <math.h>
2525 #include <stdlib.h>
26+#include <stdio.h>
2627
2728 void libtetrabz_polcmplx_1234(
2829 double g1,
@@ -268,9 +269,9 @@ void libtetrabz_polcmplx3(
268269 The latter is more accurate ?
269270 */
270271 for (i4 = 0; i4 < 4; i4++)
271- for (ir = 0; ir < 2; ir++)
272+ for (ir = 0; ir < 2; ir++) {
272273 w1[i4][ie][ir] = 0.25 / (de[i4] + e0[ie][ir]);
273-
274+ }
274275 continue;
275276
276277 for (i4 = 0; i4 < 4; i4++)
@@ -436,7 +437,7 @@ void libtetrabz_polcmplx2(
436437 */
437438 int ib, i4, j4, ir, ie, indx[4];
438439 double e[4], tsmall[4][4], v, de[4], thr, *** w2;
439-
440+
440441 w2 = (double***)malloc(4 * sizeof(double**));
441442 w2[0] = (double**)malloc(4 * ne * sizeof(double*));
442443 w2[0][0] = (double*)malloc(4 * ne * 2 * sizeof(double));
@@ -456,7 +457,7 @@ void libtetrabz_polcmplx2(
456457 for (ir = 0; ir < 2; ir++)
457458 w1[ib][i4][ie][ir] = 0.0;
458459
459- for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[i4][ib];
460+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
460461 eig_sort(4, e, indx);
461462
462463 if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
@@ -572,7 +573,7 @@ void libtetrabz_polcmplx2(
572573 }
573574 else if (e[3] <= 0.0) {
574575 for (i4 = 0; i4 < 4; i4++) {
575- de[i4] = tsmall[i4][j4] * (ej1[i4][ib] - ei1[i4]);
576+ de[i4] = ej1[ib][i4] - ei1[i4];
576577 }
577578 libtetrabz_polcmplx3(ne, e0, de, w2);
578579 for (i4 = 0; i4 < 4; i4++)
@@ -588,11 +589,11 @@ void libtetrabz_polcmplx2(
588589
589590
590591 void polcmplx(
591- int *ng,
592+ int ng[3],
592593 int nk,
593594 int nb,
594595 int ne,
595- double **bvec,
596+ double bvec[3][3],
596597 double **eig1,
597598 double **eig2,
598599 double **e0,
--- a/src/c/libtetrabz_polstat.c
+++ b/src/c/libtetrabz_polstat.c
@@ -420,7 +420,7 @@ void libtetrabz_polstat2(
420420 }
421421 else if (e[3] <= 0.0) {
422422 for (i4 = 0; i4 < 4; i4++)
423- de[i4] = ej1[i4][ib] - ei1[i4];
423+ de[i4] = ej1[ib][i4] - ei1[i4];
424424 libtetrabz_polstat3(de, w2);
425425 for (i4 = 0; i4 < 4; i4++) w1[ib][i4] += w2[i4];
426426 }
@@ -428,10 +428,10 @@ void libtetrabz_polstat2(
428428 }
429429
430430 void polstat(
431- int *ng,
431+ int ng[3],
432432 int nk,
433433 int nb,
434- double **bvec,
434+ double bvec[3][3],
435435 double **eig1,
436436 double **eig2,
437437 double ***wght
--- /dev/null
+++ b/src/c/test.c
@@ -0,0 +1,691 @@
1+/*
2+# Copyright (c) 2014 Mitsuaki Kawamura
3+#
4+# Permission is hereby granted, free of charge, to any person obtaining a
5+# copy of this software and associated documentation files (the
6+# "Software"), to deal in the Software without restriction, including
7+# without limitation the rights to use, copy, modify, merge, publish,
8+# distribute, sublicense, and/or sell copies of the Software, and to
9+# permit persons to whom the Software is furnished to do so, subject to
10+# the following conditions:
11+#
12+# The above copyright notice and this permission notice shall be included
13+# in all copies or substantial portions of the Software.
14+#
15+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
23+#include <stdio.h>
24+#include <stdlib.h>
25+#include <math.h>
26+
27+const double pi = 3.1415926535897932384626433832795;
28+void occ(int ng[3], int nk, int nb, double bvec[3][3], double** eig, double** wght);
29+void fermieng(int ng[3], int nk, int nb, double bvec[3][3], double** eig, double** wght,
30+ double nelec, int* iteration, double* ef);
31+void dos(int ng[3], int nk, int nb, int ne, double bvec[3][3], double** eig, double* e0, double*** wght);
32+void intdos(int ng[3], int nk, int nb, int ne, double bvec[3][3], double** eig, double* e0, double*** wght);
33+void dblstep(int ng[3], int nk, int nb, double bvec[3][3], double** eig1, double** eig2, double*** wght);
34+void dbldelta(int ng[3], int nk, int nb, double bvec[3][3], double** eig1, double** eig2, double*** wght);
35+void polstat(int ng[3], int nk, int nb, double bvec[3][3], double** eig1, double** eig2, double*** wght);
36+void fermigr(int ng[3], int nk, int nb, int ne, double bvec[3][3], double** eig1, double** eig2,
37+ double* e0, double**** wght);
38+void polcmplx(int ng[3], int nk, int nb, int ne, double bvec[3][3],
39+ double** eig1, double** eig2, double** e0, double***** wght);
40+
41+
42+void test_occ(
43+ int ng[3],
44+ int nk,
45+ int nb,
46+ double bvec[3][3],
47+ double vbz,
48+ double** eig,
49+ double* mat
50+) {
51+ int ib, ik;
52+ double val[2], ** wght, **eig2;
53+
54+ wght = (double**)malloc(nk * sizeof(double*));
55+ eig2 = (double**)malloc(nk * sizeof(double*));
56+ wght[0] = (double*)malloc(nk * nb * sizeof(double));
57+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
58+ for (ik = 0; ik < nk; ik++) {
59+ wght[ik] = wght[0] + ik * nb;
60+ eig2[ik] = eig2[0] + ik * nb;
61+ }
62+ for (ik = 0; ik < nk; ik++)
63+ for (ib = 0; ib < nb; ib++)
64+ eig2[ik][ib] = eig[ik][ib] - 0.5;
65+
66+ occ(ng, nk, nb, bvec, eig2, wght);
67+
68+ for (ib = 0; ib < nb; ib++) val[ib] = 0.0;
69+ for (ik = 0; ik < nk; ik++)
70+ for (ib = 0; ib < nb; ib++)
71+ val[ib] += wght[ik][ib] * mat[ik];
72+ printf("# libtetrabz_occ\n");
73+ printf(" % 15.5e % 15.5e\n", 4.0 * pi / 5.0, val[0] * vbz);
74+ printf(" % 15.5e % 15.5e\n", pi / (5.0 * sqrt(2.0)), val[1] * vbz);
75+ printf("\n");
76+
77+ free(wght[0]);
78+ free(wght);
79+ free(eig2[0]);
80+ free(eig2);
81+}
82+
83+void test_fermieng(
84+ int ng[3],
85+ int nk,
86+ int nb,
87+ double bvec[3][3],
88+ double vbz,
89+ double** eig,
90+ double* mat
91+) {
92+ int ib, ik, iteration;
93+ double val[2], ** wght, nelec, ef;
94+
95+ nelec = (4.0 * pi / 3.0 + sqrt(2.0) * pi / 3.0) / vbz;
96+
97+ wght = (double**)malloc(nk * sizeof(double*));
98+ wght[0] = (double*)malloc(nk * nb * sizeof(double));
99+ for (ik = 0; ik < nk; ik++) {
100+ wght[ik] = wght[0] + ik * nb;
101+ }
102+
103+ fermieng(ng, nk, nb, bvec, eig, wght, nelec, &iteration, &ef);
104+
105+ for (ib = 0; ib < nb; ib++) val[ib] = 0.0;
106+ for (ik = 0; ik < nk; ik++)
107+ for (ib = 0; ib < nb; ib++)
108+ val[ib] += wght[ik][ib] * mat[ik];
109+
110+ printf("# libtetrabz_fermieng\n");
111+ printf(" %15.5e %15.5e\n", 0.5, ef);
112+ printf(" %15.5e %15.5e\n", 4.0 * pi / 5.0, val[0] * vbz);
113+ printf(" %15.5e %15.5e\n", pi / (5.0 * sqrt(2.0)), val[1] * vbz);
114+ printf("\n");
115+
116+ free(wght[0]);
117+ free(wght);
118+}
119+
120+
121+void test_dos(
122+ int ng[3],
123+ int nk,
124+ int nb,
125+ double bvec[3][3],
126+ double vbz,
127+ double** eig,
128+ double* mat
129+) {
130+ int ib, ik, ne, ie;
131+ double *** wght, e0[5], val0[2][5], val[2][5];
132+
133+ ne = 5;
134+ wght = (double***)malloc(nk * sizeof(double**));
135+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
136+ wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));
137+ for (ik = 0; ik < nk; ik++) {
138+ wght[ik] = wght[0] + ik * nb;
139+ for (ib = 0; ib < nb; ib++) {
140+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
141+ }
142+ }
143+
144+ for (ie = 0; ie < ne; ie++) {
145+ e0[ie] = 0.2 * (double)(ie + 1);
146+ val0[0][ie] = 4.0 * pi * e0[ie] * e0[ie] * e0[ie];
147+ if (e0[ie] > 1.0 / sqrt(2.0))
148+ val0[1][ie] = sqrt(2.0) * pi * pow(sqrt(-1.0 + 2.0 * e0[ie] * e0[ie]), 3);
149+ else
150+ val0[1][ie] = 0.0;
151+ e0[ie] = e0[ie] * e0[ie] * 0.5;
152+ }
153+
154+ dos(ng, nk, nb, ne, bvec, eig, e0, wght);
155+
156+ for (ib = 0; ib < nb; ib++)
157+ for (ie = 0; ie < ne; ie++)
158+ val[ib][ie] = 0.0;
159+ for (ik = 0; ik < nk; ik++)
160+ for (ib = 0; ib < nb; ib++)
161+ for (ie = 0; ie < ne; ie++)
162+ val[ib][ie] += wght[ik][ib][ie] * mat[ik];
163+
164+ printf("# libtetrabz_dos\n");
165+ for (ib = 0; ib < nb; ib++)
166+ for (ie = 0; ie < ne; ie++)
167+ printf(" %15.5e %15.5e\n", val0[ib][ie], val[ib][ie] * vbz);
168+ printf("\n");
169+
170+ free(wght[0][0]);
171+ free(wght[0]);
172+ free(wght);
173+}
174+
175+void test_intdos(
176+ int ng[3],
177+ int nk,
178+ int nb,
179+ double bvec[3][3],
180+ double vbz,
181+ double** eig,
182+ double* mat
183+) {
184+ int ib, ik, ne, ie;
185+ double*** wght, e0[5], val0[2][5], val[2][5];
186+
187+ ne = 5;
188+ wght = (double***)malloc(nk * sizeof(double**));
189+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
190+ wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));
191+ for (ik = 0; ik < nk; ik++) {
192+ wght[ik] = wght[0] + ik * nb;
193+ for (ib = 0; ib < nb; ib++) {
194+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
195+ }
196+ }
197+
198+ for (ie = 0; ie < ne; ie++) {
199+ e0[ie] = 0.2 * (double)(ie + 1);
200+ val0[0][ie] = 4.0 * pi * pow(e0[ie], 5) / 5.0;
201+ if (e0[ie] > 1.0 / sqrt(2.0))
202+ val0[1][ie] = pi * pow(sqrt(-1.0 + 2.0 * pow(e0[ie], 2)), 5) / (5.0 * sqrt(2.0));
203+ else
204+ val0[1][ie] = 0.0;
205+ e0[ie] = pow(e0[ie], 2) * 0.5;
206+ }
207+
208+ intdos(ng, nk, nb, ne, bvec, eig, e0, wght);
209+
210+ for (ib = 0; ib < nb; ib++)
211+ for (ie = 0; ie < ne; ie++)
212+ val[ib][ie] = 0.0;
213+ for (ik = 0; ik < nk; ik++)
214+ for (ib = 0; ib < nb; ib++)
215+ for (ie = 0; ie < ne; ie++)
216+ val[ib][ie] += wght[ik][ib][ie] * mat[ik];
217+
218+ printf("# libtetrabz_intdos\n");
219+ for (ib = 0; ib < nb; ib++)
220+ for (ie = 0; ie < ne; ie++)
221+ printf(" %15.5e %15.5e\n", val0[ib][ie], val[ib][ie] * vbz);
222+ printf("\n");
223+
224+ free(wght[0][0]);
225+ free(wght[0]);
226+ free(wght);
227+}
228+
229+
230+void test_dblstep(
231+ int ng[3],
232+ int nk,
233+ int nb,
234+ double bvec[3][3],
235+ double vbz,
236+ double** eig1,
237+ double** eig2,
238+ double* mat
239+) {
240+ int ib, ik, jb;
241+ double*** wght, val[2][2], ** eig12, ** eig22;
242+
243+ wght = (double***)malloc(nk * sizeof(double**));
244+ eig12 = (double**)malloc(nk * sizeof(double*));
245+ eig22 = (double**)malloc(nk * sizeof(double*));
246+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
247+ eig12[0] = (double*)malloc(nk * nb * sizeof(double));
248+ eig22[0] = (double*)malloc(nk * nb * sizeof(double));
249+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
250+ for (ik = 0; ik < nk; ik++) {
251+ wght[ik] = wght[0] + ik * nb;
252+ eig12[ik] = eig12[0] + ik * nb;
253+ eig22[ik] = eig22[0] + ik * nb;
254+ for (ib = 0; ib < nb; ib++) {
255+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
256+ }
257+ }
258+ for (ik = 0; ik < nk; ik++)
259+ for (ib = 0; ib < nb; ib++) {
260+ eig12[ik][ib] = eig1[ik][ib] - 0.5;
261+ eig22[ik][ib] = eig2[ik][ib] - 0.5;
262+ }
263+
264+ dblstep(ng, nk, nb, bvec, eig12, eig22, wght);
265+
266+ for (ib = 0; ib < nb; ib++)
267+ for (jb = 0; jb < nb; jb++)
268+ val[ib][jb] = 0.0;
269+ for (ik = 0; ik < nk; ik++)
270+ for (ib = 0; ib < nb; ib++)
271+ for (jb = 0; jb < nb; jb++)
272+ val[ib][jb] += wght[ik][ib][jb] * mat[ik];
273+
274+ printf("# libtetrabz_dblstep\n");
275+ printf(" %15.5e %15.5e\n", 49.0 * pi / 320.0, val[0][0] * vbz);
276+ printf(" %15.5e %15.5e\n", 0.0, val[0][1] * vbz);
277+ printf(" %15.5e %15.5e\n", pi * (512.0 * sqrt(2.0) - 319.0) / 10240.0, val[1][0] * vbz);
278+ printf(" %15.5e %15.5e\n", 0.0, val[1][1] * vbz);
279+ printf("\n");
280+
281+ free(wght[0][0]);
282+ free(wght[0]);
283+ free(wght);
284+ free(eig12[0]);
285+ free(eig12);
286+ free(eig22[0]);
287+ free(eig22);
288+}
289+
290+
291+void test_dbldelta(
292+ int ng[3],
293+ int nk,
294+ int nb,
295+ double bvec[3][3],
296+ double vbz,
297+ double** eig1,
298+ double** eig2,
299+ double* mat
300+) {
301+ int ib, ik, jb;
302+ double*** wght, val[2][2], ** eig12, ** eig22;
303+
304+ wght = (double***)malloc(nk * sizeof(double**));
305+ eig12 = (double**)malloc(nk * sizeof(double*));
306+ eig22 = (double**)malloc(nk * sizeof(double*));
307+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
308+ eig12[0] = (double*)malloc(nk * nb * sizeof(double));
309+ eig22[0] = (double*)malloc(nk * nb * sizeof(double));
310+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
311+ for (ik = 0; ik < nk; ik++) {
312+ wght[ik] = wght[0] + ik * nb;
313+ eig12[ik] = eig12[0] + ik * nb;
314+ eig22[ik] = eig22[0] + ik * nb;
315+ for (ib = 0; ib < nb; ib++) {
316+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
317+ }
318+ }
319+ for (ik = 0; ik < nk; ik++)
320+ for (ib = 0; ib < nb; ib++) {
321+ eig12[ik][ib] = eig1[ik][ib] - 0.5;
322+ eig22[ik][ib] = eig2[ik][ib] - 0.5;
323+ }
324+
325+ dbldelta(ng, nk, nb, bvec, eig12, eig22, wght);
326+
327+ for (ib = 0; ib < nb; ib++)
328+ for (jb = 0; jb < nb; jb++)
329+ val[ib][jb] = 0.0;
330+ for (ik = 0; ik < nk; ik++)
331+ for (ib = 0; ib < nb; ib++)
332+ for (jb = 0; jb < nb; jb++)
333+ val[ib][jb] += wght[ik][ib][jb] * mat[ik];
334+ printf("# libtetrabz_dbldelta\n");
335+ printf(" %15.5e %15.5e\n", 2.0 * pi, val[0][0] * vbz);
336+ printf(" %15.5e %15.5e\n", 0.0, val[0][1] * vbz);
337+ printf(" %15.5e %15.5e\n", pi, val[1][0] * vbz);
338+ printf(" %15.5e %15.5e\n", 0.0, val[1][1] * vbz);
339+ printf("\n");
340+
341+ free(wght[0][0]);
342+ free(wght[0]);
343+ free(wght);
344+ free(eig12[0]);
345+ free(eig12);
346+ free(eig22[0]);
347+ free(eig22);
348+}
349+
350+
351+void test_polstat(
352+ int ng[3],
353+ int nk,
354+ int nb,
355+ double bvec[3][3],
356+ double vbz,
357+ double** eig1,
358+ double** eig2,
359+ double* mat
360+) {
361+ int ib, ik, jb;
362+ double*** wght, val[2][2], ** eig12, ** eig22;
363+
364+ wght = (double***)malloc(nk * sizeof(double**));
365+ eig12 = (double**)malloc(nk * sizeof(double*));
366+ eig22 = (double**)malloc(nk * sizeof(double*));
367+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
368+ eig12[0] = (double*)malloc(nk * nb * sizeof(double));
369+ eig22[0] = (double*)malloc(nk * nb * sizeof(double));
370+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
371+ for (ik = 0; ik < nk; ik++) {
372+ wght[ik] = wght[0] + ik * nb;
373+ eig12[ik] = eig12[0] + ik * nb;
374+ eig22[ik] = eig22[0] + ik * nb;
375+ for (ib = 0; ib < nb; ib++) {
376+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
377+ }
378+ }
379+ for (ik = 0; ik < nk; ik++)
380+ for (ib = 0; ib < nb; ib++) {
381+ eig12[ik][ib] = eig1[ik][ib] - 0.5;
382+ eig22[ik][ib] = eig2[ik][ib] - 0.5;
383+ }
384+
385+ polstat(ng, nk, nb, bvec, eig12, eig22, wght);
386+
387+ for (ib = 0; ib < nb; ib++)
388+ for (jb = 0; jb < nb; jb++)
389+ val[ib][jb] = 0.0;
390+ for (ik = 0; ik < nk; ik++)
391+ for (ib = 0; ib < nb; ib++)
392+ for (jb = 0; jb < nb; jb++)
393+ val[ib][jb] += wght[ik][ib][jb] * mat[ik];
394+ printf("# libtetrabz_polstat");
395+ printf(" %15.5e %15.5e", pi * (68.0 + 45.0 * log(3.0)) / 96.0, val[0][0] * vbz);
396+ printf(" %15.5e %15.5e", pi * 8.0 / 5.0, val[0][1] * vbz);
397+ printf(" %15.5e %15.5e", pi * (228.0 + 22.0 * sqrt(2.0) - 96.0 * log(2.0)
398+ + 192.0 * log(4.0 + sqrt(2.0))
399+ - 3.0 * log(1.0 + 2.0 * sqrt(2.0))) / 1536.0,
400+ val[1][0] * vbz);
401+ printf(" %15.5e %15.5e", pi * sqrt(8.0) / 5.0, val[1][1] * vbz);
402+ printf("\n");
403+
404+ free(wght[0][0]);
405+ free(wght[0]);
406+ free(wght);
407+ free(eig12[0]);
408+ free(eig12);
409+ free(eig22[0]);
410+ free(eig22);
411+}
412+
413+
414+void test_fermigr(
415+ int ng[3],
416+ int nk,
417+ int nb,
418+ double bvec[3][3],
419+ double vbz,
420+ double** eig1,
421+ double** eig2,
422+ double* mat
423+) {
424+ int ib, ik, jb, ne, ie;
425+ double**** wght, val[2][2][3], val0[2][2][3], ** eig12, ** eig22, e0[3];
426+
427+ ne = 3;
428+
429+ wght = (double****)malloc(nk * sizeof(double***));
430+ eig12 = (double**)malloc(nk * sizeof(double*));
431+ eig22 = (double**)malloc(nk * sizeof(double*));
432+ wght[0] = (double***)malloc(nk * nb * sizeof(double**));
433+ eig12[0] = (double*)malloc(nk * nb * sizeof(double));
434+ eig22[0] = (double*)malloc(nk * nb * sizeof(double));
435+ wght[0][0] = (double**)malloc(nk * nb * nb * sizeof(double*));
436+ wght[0][0][0] = (double*)malloc(nk * nb * nb * ne * sizeof(double));
437+ for (ik = 0; ik < nk; ik++) {
438+ wght[ik] = wght[0] + ik * nb;
439+ eig12[ik] = eig12[0] + ik * nb;
440+ eig22[ik] = eig22[0] + ik * nb;
441+ for (ib = 0; ib < nb; ib++) {
442+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
443+ for (jb = 0; jb < nb; jb++) {
444+ wght[ik][ib][jb] = wght[0][0][0] + ik * nb * nb * ne + ib * nb * ne + jb * ne;
445+ }
446+ }
447+ }
448+ for (ik = 0; ik < nk; ik++)
449+ for (ib = 0; ib < nb; ib++) {
450+ eig12[ik][ib] = eig1[ik][ib] - 0.5;
451+ eig22[ik][ib] = eig2[ik][ib] - 0.5;
452+ }
453+
454+ e0[0] = 1.0 / 3.0;
455+ e0[1] = 2.0 / 3.0;
456+ e0[2] = 1.0;
457+ val0[0][0][0] = 4.0 * pi / 9.0;
458+ val0[0][0][1] = 1295.0 * pi / 2592.0;
459+ val0[0][0][2] = 15.0 * pi / 32.0;
460+ val0[1][0][0] = 5183.0 * pi / 41472.0;
461+ val0[1][0][1] = 4559.0 * pi / 41472.0;
462+ val0[1][0][2] = 0.0;
463+ for (ib = 0; ib < nb; ib++)
464+ for (ie = 0; ie < ne; ie++)
465+ val0[ib][1][ie] = 0.0;
466+
467+ fermigr(ng, nk, nb, ne, bvec, eig12, eig22, e0, wght);
468+
469+ for (ib = 0; ib < nb; ib++)
470+ for (jb = 0; jb < nb; jb++)
471+ for (ie = 0; ie < ne; ie++)
472+ val[ib][jb][ie] = 0.0;
473+ for (ik = 0; ik < nk; ik++)
474+ for (ib = 0; ib < nb; ib++)
475+ for (jb = 0; jb < nb; jb++)
476+ for (ie = 0; ie < ne; ie++)
477+ val[ib][jb][ie] += wght[ik][ib][jb][ie] * mat[ik];
478+
479+ printf("# libtetrabz_fermigr\n");
480+ for (ib = 0; ib < nb; ib++)
481+ for (jb = 0; jb < nb; jb++)
482+ for (ie = 0; ie < ne; ie++)
483+ printf(" %15.5e %15.5e\n", val0[ib][jb][ie], val[ib][jb][ie] * vbz);
484+ printf("\n");
485+
486+ free(wght[0][0][0]);
487+ free(wght[0][0]);
488+ free(wght[0]);
489+ free(wght);
490+ free(eig12[0]);
491+ free(eig12);
492+ free(eig22[0]);
493+ free(eig22);
494+}
495+
496+void test_polcmplx(
497+ int ng[3],
498+ int nk,
499+ int nb,
500+ double bvec[3][3],
501+ double vbz,
502+ double** eig1,
503+ double** eig2,
504+ double* mat
505+) {
506+ int ib, ik, jb, ne, ie, ir;
507+ double***** wght, val[2][2][3][2], val0[2][2][3][2], ** eig12, ** eig22, **e0;
508+
509+ ne = 3;
510+
511+ e0 = (double**)malloc(ne * sizeof(double*));
512+ e0[0] = (double*)malloc(ne * 2 * sizeof(double));
513+ for (ie = 0; ie < ne; ie++)
514+ e0[ie] = e0[0] + ie * 2;
515+
516+ wght = (double*****)malloc(nk * sizeof(double****));
517+ eig12 = (double**)malloc(nk * sizeof(double*));
518+ eig22 = (double**)malloc(nk * sizeof(double*));
519+ wght[0] = (double****)malloc(nk * nb * sizeof(double***));
520+ eig12[0] = (double*)malloc(nk * nb * sizeof(double));
521+ eig22[0] = (double*)malloc(nk * nb * sizeof(double));
522+ wght[0][0] = (double***)malloc(nk * nb * nb * sizeof(double**));
523+ wght[0][0][0] = (double**)malloc(nk * nb * nb * ne * sizeof(double*));
524+ wght[0][0][0][0] = (double*)malloc(nk * nb * nb * ne * 2 * sizeof(double));
525+ for (ik = 0; ik < nk; ik++) {
526+ wght[ik] = wght[0] + ik * nb;
527+ eig12[ik] = eig12[0] + ik * nb;
528+ eig22[ik] = eig22[0] + ik * nb;
529+ for (ib = 0; ib < nb; ib++) {
530+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
531+ for (jb = 0; jb < nb; jb++) {
532+ wght[ik][ib][jb] = wght[0][0][0] + ik * nb * nb * ne + ib * nb * ne + jb * ne;
533+ for (ie = 0; ie < ne; ie++) {
534+ wght[ik][ib][jb][ie] = wght[0][0][0][0] + ik * nb * nb * ne * 2
535+ + ib * nb * ne * 2 + jb * ne * 2 + ie * 2;
536+ }
537+ }
538+ }
539+ }
540+ for (ik = 0; ik < nk; ik++)
541+ for (ib = 0; ib < nb; ib++) {
542+ eig12[ik][ib] = eig1[ik][ib] - 0.5;
543+ eig22[ik][ib] = eig2[ik][ib] - 0.5;
544+ }
545+
546+ e0[0][0] = -2.0;
547+ e0[0][1] = 1.0;
548+ e0[1][0] = 0.0;
549+ e0[1][1] = 2.0;
550+ e0[2][0] = 1.0;
551+ e0[2][1] = -0.5;
552+ val0[0][0][0][0] = -0.838243341280338;
553+ val0[0][0][0][1] = -0.734201894333234;
554+ val0[0][0][1][0] = 0.270393588876530;
555+ val0[0][0][1][1] = -0.771908416949610;
556+ val0[0][0][2][0] = 0.970996830573510;
557+ val0[0][0][2][1] = 0.302792326476720;
558+ val0[1][0][0][0] = -0.130765724778920;
559+ val0[1][0][0][1] = -0.087431218706638;
560+ val0[1][0][1][0] = 0.030121954547245;
561+ val0[1][0][1][1] = -0.135354254293510;
562+ val0[1][0][2][0] = 0.178882244951203;
563+ val0[1][0][2][1] = 0.064232167683425;
564+
565+ for (ie = 0; ie < ne; ie++) {
566+ for (ir = 0; ir < 2; ir++) {
567+ val0[0][1][ie][ir] = (8.0 * pi) / (5.0 * (1.0 + 2.0 * e0[ie][ir]));
568+ val0[1][1][ie][ir] = (sqrt(8.0) * pi) / (5.0 * (1.0 + 4.0 * e0[ie][ir]));
569+ }
570+ }
571+
572+ polcmplx(ng, nk, nb, ne, bvec, eig12, eig22, &e0[0], wght);
573+
574+ for (ib = 0; ib < nb; ib++)
575+ for (jb = 0; jb < nb; jb++)
576+ for (ie = 0; ie < ne; ie++)
577+ for (ir = 0; ir < 2; ir++)
578+ val[ib][jb][ie][ir] = 0.0;
579+ for (ik = 0; ik < nk; ik++)
580+ for (ib = 0; ib < nb; ib++)
581+ for (jb = 0; jb < nb; jb++)
582+ for (ie = 0; ie < ne; ie++)
583+ for (ir = 0; ir < 2; ir++)
584+ val[ib][jb][ie][ir] += wght[ik][ib][jb][ie][ir] * mat[ik];
585+
586+
587+ printf("# libtetrabz_polcmplx\n");
588+ for (ib = 0; ib < nb; ib++)
589+ for (jb = 0; jb < nb; jb++)
590+ for (ie = 0; ie < ne; ie++) {
591+ printf(" %15.5e %15.5e\n", val0[ib][jb][ie][0], val[ib][jb][ie][0] * vbz);
592+ printf(" %15.5e %15.5e\n", val0[ib][jb][ie][1], val[ib][jb][ie][1] * vbz);
593+ }
594+ printf("\n");
595+ free(wght[0][0][0]);
596+ free(wght[0][0]);
597+ free(wght[0]);
598+ free(wght);
599+ free(eig12[0]);
600+ free(eig12);
601+ free(eig22[0]);
602+ free(eig22);
603+ free(e0[0]);
604+ free(e0);
605+}
606+
607+int main()
608+{
609+ int ng0, ng[3], nb, i3, j3, nk, ik, i0, i1, i2;
610+ double bvec[3][3], vbz, ** eig1, ** eig2, * mat, kvec[3];
611+
612+ ng0 = 8;
613+ for (i3 = 0; i3 < 3; i3++) ng[i3] = ng0;
614+ nk = ng[0] * ng[1] * ng[2];
615+ nb = 2;
616+ for (i3 = 0; i3 < 3; i3++)
617+ for (j3 = 0; j3 < 3; j3++)
618+ bvec[i3][j3] = 0.0;
619+ for (i3 = 0; i3 < 3; i3++) bvec[i3][i3] = 3.0;
620+ vbz = 27.0;/*abs(numpy.linalg.det(bvec));*/
621+
622+ mat = (double*)malloc(nk * sizeof(double));
623+ eig1 = (double**)malloc(nk * sizeof(double*));
624+ eig2 = (double**)malloc(nk * sizeof(double*));
625+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
626+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
627+ for (ik = 0; ik < nk; ik++) {
628+ eig1[ik] = eig1[0] + ik * nb;
629+ eig2[ik] = eig2[0] + ik * nb;
630+ }
631+
632+ for (i0 = 0; i0 < ng[0]; i0++) {
633+ for (i1 = 0; i1 < ng[1]; i1++) {
634+ for (i2 = 0; i2 < ng[2]; i2++) {
635+
636+ kvec[0] = (double)i0 / (double)ng[0];
637+ kvec[1] = (double)i1 / (double)ng[1];
638+ kvec[2] = (double)i2 / (double)ng[2];
639+ for (i3 = 0; i3 < 3; i3++)
640+ if (kvec[i3] > 0.5) kvec[i3] += -1.0;
641+ for (i3 = 0; i3 < 3; i3++) kvec[i3] = kvec[i3] * bvec[i3][i3];/*kvec.dot(bvec)*/
642+
643+ ik = i2 + i1 * ng[2] + i0 * ng[1] * ng[2];
644+ eig1[ik][0] = 0.5 * (kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2]);
645+ eig1[ik][1] = eig1[ik][0] + 0.25;
646+
647+ kvec[0] += 1.0;
648+ eig2[ik][0] = 0.5 * (kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2]);
649+ eig2[ik][1] = eig1[ik][0] + 0.5;
650+ }
651+ }
652+ }
653+
654+ for (i0 = 0; i0 < ng[0]; i0++) {
655+ for (i1 = 0; i1 < ng[1]; i1++) {
656+ for (i2 = 0; i2 < ng[2]; i2++) {
657+
658+ kvec[0] = (double)i0 / (double)ng[0];
659+ kvec[1] = (double)i1 / (double)ng[1];
660+ kvec[2] = (double)i2 / (double)ng[2];
661+ for (i3 = 0; i3 < 3; i3++)
662+ if (kvec[i3] > 0.5) kvec[i3] += -1.0;
663+ for (i3 = 0; i3 < 3; i3++) kvec[i3] = kvec[i3] * bvec[i3][i3];/*kvec.dot(bvec)*/
664+
665+ ik = i2 + i1 * ng[2] + i0 * ng[1] * ng[2];
666+ mat[ik] = kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2];
667+ }
668+ }
669+ }
670+
671+ printf("# Ideal Result\n");
672+ /*
673+ test_occ(ng, nk, nb, bvec, vbz, eig1, mat);
674+
675+ test_fermieng(ng, nk, nb, bvec, vbz, eig1, mat);
676+
677+ test_dos(ng, nk, nb, bvec, vbz, eig1, mat);
678+
679+ test_intdos(ng, nk, nb, bvec, vbz, eig1, mat);
680+ */
681+ test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
682+ /*
683+ test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
684+
685+ test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
686+
687+ test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
688+
689+ test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
690+ */
691+}
Show on old repository browser