• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revisionfede19916576a39166f82368081f4886ed547d24 (tree)
Zeit2022-03-24 01:25:06
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/src/c/libtetrabz.c
@@ -0,0 +1,4116 @@
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+#define PY_SSIZE_T_CLEAN
24+#include <Python.h>
25+#include <stdio.h>
26+#include <stdlib.h>
27+#include <math.h>
28+
29+void libtetrabz_initialize(
30+ int ng[3],
31+ double bvec[3][3],
32+ double wlsm[20][4],
33+ int **ikv
34+ )
35+ {
36+ /*
37+ define shortest diagonal line & define type of tetragonal
38+ */
39+ int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
40+ double bvec2[3][3], bvec3[4][3], bnorm[4];
41+ double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
42+ {0.0, 1440.0, 0.0, 30.0},
43+ {30.0, 0.0, 1440.0, 0.0},
44+ {0.0, 30.0, 0.0, 1440.0} },
45+ wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
46+ {-28.0, -38.0, 7.0, 17.0},
47+ {17.0, -28.0, -38.0, 7.0},
48+ {7.0, 17.0, -28.0, -38.0}},
49+ wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
50+ {9.0, -56.0, 9.0, -46.0},
51+ {-46.0, 9.0, -56.0, 9.0},
52+ {9.0, -46.0, 9.0, -56.0}},
53+ wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
54+ {7.0, -38.0, -28.0, 17.0},
55+ {17.0, 7.0, -38.0, -28.0},
56+ {-28.0, 17.0, 7.0, -38.0}},
57+ wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
58+ {-18.0, -18.0, -18.0, 12.0},
59+ {12.0, -18.0, -18.0, -18.0},
60+ {-18.0, 12.0, -18.0, -18.0}};
61+
62+ for (i1 = 0; i1 < 3; i1++)
63+ for (i2 = 0; i2 < 3; i2++)
64+ bvec2[i1][i2] = bvec[i1][i2] / (double) ng[i1];
65+
66+ for (i1 = 0; i1 < 3; i1++) {
67+ bvec3[0][i1] = -bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
68+ bvec3[1][i1] = bvec2[0][i1] - bvec2[1][i1] + bvec2[2][i1];
69+ bvec3[2][i1] = bvec2[0][i1] + bvec2[1][i1] - bvec2[2][i1];
70+ bvec3[3][i1] = bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
71+ }
72+ /*
73+ length of delta bvec
74+ */
75+ for (i1 = 0; i1 < 4; i1++) {
76+ bnorm[i1] = 0.0;
77+ for (i2 = 0; i2 < 3; i2++)
78+ bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
79+ }
80+ itype = 0;
81+ for (i1 = 1; i1 < 4; i1++)
82+ if (bnorm[i1] < bnorm[itype]) itype = i1;
83+ /*
84+ start & last
85+ */
86+ for (i0 = 0; i0 < 4; i0++) {
87+ ivvec0[i0] = 0;
88+ for (i1 = 0; i1 < 4; i1++)divvec[i0][i1] = 0;
89+ divvec[i0][i0] = 1;
90+ }
91+ ivvec0[itype] = 1;
92+ divvec[itype][itype] = -1;
93+ /*
94+ Corners of tetrahedron
95+ */
96+ it = 0;
97+ for (i0 = 0; i0 < 3; i0++) {
98+ for (i1 = 0; i1 < 3; i1++) {
99+ if (i1 == i0) continue;
100+ for (i2 = 0; i2 < 3; i2++) {
101+ if (i2 == i1 || i2 == i0) continue;
102+
103+ for (i3 = 0; i3 < 3; i3++) {
104+ ivvec[it][0][i3] = ivvec0[i3];
105+ ivvec[it][1][i3] = ivvec[it][0][i3] + divvec[i0][i3];
106+ ivvec[it][2][i3] = ivvec[it][1][i3] + divvec[i1][i3];
107+ ivvec[it][3][i3] = ivvec[it][2][i3] + divvec[i2][i3];
108+ }
109+
110+ it += 1;
111+ }
112+ }
113+ }
114+ /*
115+ Additional points
116+ */
117+ for (i1 = 0; i1 < 6; i1++) {
118+ for (i2 = 0; i2 < 3; i2++) {
119+ ivvec[i1][4][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][1][i2];
120+ ivvec[i1][5][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][2][i2];
121+ ivvec[i1][6][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][3][i2];
122+ ivvec[i1][7][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][0][i2];
123+
124+ ivvec[i1][8][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][2][i2];
125+ ivvec[i1][9][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][3][i2];
126+ ivvec[i1][10][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][0][i2];
127+ ivvec[i1][11][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][1][i2];
128+
129+ ivvec[i1][12][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][3][i2];
130+ ivvec[i1][13][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][0][i2];
131+ ivvec[i1][14][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][1][i2];
132+ ivvec[i1][15][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][2][i2];
133+
134+ ivvec[i1][16][i2] = ivvec[i1][3][i2] - ivvec[i1][0][i2] + ivvec[i1][1][i2];
135+ ivvec[i1][17][i2] = ivvec[i1][0][i2] - ivvec[i1][1][i2] + ivvec[i1][2][i2];
136+ ivvec[i1][18][i2] = ivvec[i1][1][i2] - ivvec[i1][2][i2] + ivvec[i1][3][i2];
137+ ivvec[i1][19][i2] = ivvec[i1][2][i2] - ivvec[i1][3][i2] + ivvec[i1][0][i2];
138+ }
139+ }
140+
141+ for (i1 = 0; i1 < 4; i1++) {
142+ for (i2 = 0; i2 < 4; i2++) {
143+ wlsm[i2][i1] = wlsm1[i1][i2] /= 1260.0;
144+ wlsm[i2 + 4][i1] = wlsm2[i1][i2] /= 1260.0;
145+ wlsm[i2 + 8][i1] = wlsm3[i1][i2] /= 1260.0;
146+ wlsm[i2 + 12][i1] = wlsm4[i1][i2] /= 1260.0;
147+ wlsm[i2 + 16][i1] = wlsm5[i1][i2] /= 1260.0;
148+ }
149+ }
150+ /*
151+ k-index for energy
152+ */
153+ nt = 0;
154+ for (i2 = 0; i2 < ng[2]; i2++) {
155+ for (i1 = 0; i1 < ng[1]; i1++) {
156+ for (i0 = 0; i0 < ng[0]; i0++) {
157+
158+ for (it = 0; it < 6; it++) {
159+
160+ for (ii = 0; ii < 20; ii++) {
161+ ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
162+ ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
163+ ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
164+ for (i3 = 0; i3 < 3; i3++) if (ikv0[i3] < 0) ikv0[i3] += ng[i3];
165+ ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
166+ }
167+ nt += 1;
168+ }
169+ }
170+ }
171+ }
172+}
173+/*
174+Cut small tetrahedron A1
175+*/
176+void libtetrabz_tsmall_a1(
177+ double *e,
178+ double e0,
179+ double *v,
180+ double tsmall[4][4]
181+) {
182+ double a10, a20, a30;
183+ a10 = (e0 - e[0]) / (e[1] - e[0]);
184+ a20 = (e0 - e[0]) / (e[2] - e[0]);
185+ a30 = (e0 - e[0]) / (e[3] - e[0]);
186+
187+ *v = a10 * a20 * a30;
188+
189+ tsmall[0][0] = 1.0;
190+ tsmall[0][1] = 1.0 - a10;
191+ tsmall[0][2] = 1.0 - a20;
192+ tsmall[0][3] = 1.0 - a30;
193+ tsmall[1][0] = 0.0;
194+ tsmall[1][1] = a10;
195+ tsmall[1][2] = 0.0;
196+ tsmall[1][3] = 0.0;
197+ tsmall[2][0] = 0.0;
198+ tsmall[2][1] = 0.0;
199+ tsmall[2][2] = a20;
200+ tsmall[2][3] = 0.0;
201+ tsmall[3][0] = 0.0;
202+ tsmall[3][1] = 0.0;
203+ tsmall[3][2] = 0.0;
204+ tsmall[3][3] = a30;
205+}
206+/*
207+Cut small tetrahedron B1
208+*/
209+void libtetrabz_tsmall_b1(
210+ double *e,
211+ double e0,
212+ double *v,
213+ double tsmall[4][4]
214+)
215+{
216+ double a13, a20, a30;
217+ a13 = (e0 - e[3]) / (e[1] - e[3]);
218+ a20 = (e0 - e[0]) / (e[2] - e[0]);
219+ a30 = (e0 - e[0]) / (e[3] - e[0]);
220+
221+ *v = a20 * a30 * a13;
222+
223+ tsmall[0][0] = 1.0;
224+ tsmall[0][1] = 1.0 - a20;
225+ tsmall[0][2] = 1.0 - a30;
226+ tsmall[0][3] = 0.0;
227+ tsmall[1][0] = 0.0;
228+ tsmall[1][1] = 0.0;
229+ tsmall[1][2] = 0.0;
230+ tsmall[1][3] = a13;
231+ tsmall[2][0] = 0.0;
232+ tsmall[2][1] = a20;
233+ tsmall[2][2] = 0.0;
234+ tsmall[2][3] = 0.0;
235+ tsmall[3][0] = 0.0;
236+ tsmall[3][1] = 0.0;
237+ tsmall[3][2] = a30;
238+ tsmall[3][3] = 1.0 - a13;
239+}
240+/*
241+Cut small tetrahedron B2
242+*/
243+void libtetrabz_tsmall_b2(
244+ double *e,
245+ double e0,
246+ double *v,
247+ double tsmall[4][4]
248+)
249+{
250+ double a21, a31;
251+ a21 = (e0 - e[1]) / (e[2] - e[1]);
252+ a31 = (e0 - e[1]) / (e[3] - e[1]);
253+
254+ *v = a21 * a31;
255+
256+ tsmall[0][0] = 1.0;
257+ tsmall[0][1] = 0.0;
258+ tsmall[0][2] = 0.0;
259+ tsmall[0][3] = 0.0;
260+ tsmall[1][0] = 0.0;
261+ tsmall[1][1] = 1.0;
262+ tsmall[1][2] = 1.0 - a21;
263+ tsmall[1][3] = 1.0 - a31;
264+ tsmall[2][0] = 0.0;
265+ tsmall[2][1] = 0.0;
266+ tsmall[2][2] = a21;
267+ tsmall[2][3] = 0.0;
268+ tsmall[3][0] = 0.0;
269+ tsmall[3][1] = 0.0;
270+ tsmall[3][2] = 0.0;
271+ tsmall[3][3] = a31;
272+}
273+/*
274+Cut small tetrahedron B3
275+*/
276+void libtetrabz_tsmall_b3(
277+ double *e,
278+ double e0,
279+ double *v,
280+ double tsmall[4][4]
281+)
282+{
283+ double a12, a20, a31;
284+ a12 = (e0 - e[2]) / (e[1] - e[2]);
285+ a20 = (e0 - e[0]) / (e[2] - e[0]);
286+ a31 = (e0 - e[1]) / (e[3] - e[1]);
287+
288+ *v = a12 * a20 * a31;
289+
290+ tsmall[0][0] = 1.0;
291+ tsmall[0][1] = 1.0 - a20;
292+ tsmall[0][2] = 0.0;
293+ tsmall[0][3] = 0.0;
294+ tsmall[1][0] = 0.0;
295+ tsmall[1][1] = 0.0;
296+ tsmall[1][2] = a12;
297+ tsmall[1][3] = 1.0 - a31;
298+ tsmall[2][0] = 0.0;
299+ tsmall[2][1] = a20;
300+ tsmall[2][2] = 1.0 - a12;
301+ tsmall[2][3] = 0.0;
302+ tsmall[3][0] = 0.0;
303+ tsmall[3][1] = 0.0;
304+ tsmall[3][2] = 0.0;
305+ tsmall[3][3] = a31;
306+}
307+/*
308+Cut small tetrahedron C1
309+*/
310+void libtetrabz_tsmall_c1(
311+ double *e,
312+ double e0,
313+ double *v,
314+ double tsmall[4][4]
315+)
316+{
317+ double a32;
318+ a32 = (e0 - e[2]) / (e[3] - e[2]);
319+
320+ *v = a32;
321+
322+ tsmall[0][0] = 1.0;
323+ tsmall[0][1] = 0.0;
324+ tsmall[0][2] = 0.0;
325+ tsmall[0][3] = 0.0;
326+ tsmall[1][0] = 0.0;
327+ tsmall[1][1] = 1.0;
328+ tsmall[1][2] = 0.0;
329+ tsmall[1][3] = 0.0;
330+ tsmall[2][0] = 0.0;
331+ tsmall[2][1] = 0.0;
332+ tsmall[2][2] = 1.0;
333+ tsmall[2][3] = 1.0 - a32;
334+ tsmall[3][0] = 0.0;
335+ tsmall[3][1] = 0.0;
336+ tsmall[3][2] = 0.0;
337+ tsmall[3][3] = a32;
338+}
339+/*
340+Cut small tetrahedron C2
341+*/
342+void libtetrabz_tsmall_c2(
343+ double *e,
344+ double e0,
345+ double *v,
346+ double tsmall[4][4]
347+)
348+{
349+ double a23, a31;
350+ a23 = (e0 - e[3]) / (e[2] - e[3]);
351+ a31 = (e0 - e[1]) / (e[3] - e[1]);
352+
353+ *v = a23 * a31;
354+
355+ tsmall[0][0] = 1.0;
356+ tsmall[0][1] = 0.0;
357+ tsmall[0][2] = 0.0;
358+ tsmall[0][3] = 0.0;
359+ tsmall[1][0] = 0.0;
360+ tsmall[1][1] = 1.0;
361+ tsmall[1][2] = 1.0 - a31;
362+ tsmall[1][3] = 0.0;
363+ tsmall[2][0] = 0.0;
364+ tsmall[2][1] = 0.0;
365+ tsmall[2][2] = 0.0;
366+ tsmall[2][3] = a23;
367+ tsmall[3][0] = 0.0;
368+ tsmall[3][1] = 0.0;
369+ tsmall[3][2] = a31;
370+ tsmall[3][3] = 1.0 - a23;
371+}
372+/*
373+Cut small tetrahedron C3
374+*/
375+void libtetrabz_tsmall_c3(
376+ double *e,
377+ double e0,
378+ double *v,
379+ double tsmall[4][4]
380+)
381+{
382+ double a23, a13, a30;
383+ a23 = (e0 - e[3]) / (e[2] - e[3]);
384+ a13 = (e0 - e[3]) / (e[1] - e[3]);
385+ a30 = (e0 - e[0]) / (e[3] - e[0]);
386+
387+ *v = a23 * a13 * a30;
388+
389+ tsmall[0][0] = 1.0;
390+ tsmall[0][1] = 1.0 - a30;
391+ tsmall[0][2] = 0.0;
392+ tsmall[0][3] = 0.0;
393+ tsmall[1][0] = 0.0;
394+ tsmall[1][1] = 0.0;
395+ tsmall[1][2] = a13;
396+ tsmall[1][3] = 0.0;
397+ tsmall[2][0] = 0.0;
398+ tsmall[2][1] = 0.0;
399+ tsmall[2][2] = 0.0;
400+ tsmall[2][3] = a23;
401+ tsmall[3][0] = 0.0;
402+ tsmall[3][1] = a30;
403+ tsmall[3][2] = 1.0 - a13;
404+ tsmall[3][3] = 1.0 - a23;
405+}
406+/*
407+Cut triangle A1
408+*/
409+void libtetrabz_triangle_a1(
410+ double *e,
411+ double e0,
412+ double *v,
413+ double tsmall[4][3]
414+)
415+{
416+ double a10, a20, a30;
417+ a10 = (e0 - e[0]) / (e[1] - e[0]);
418+ a20 = (e0 - e[0]) / (e[2] - e[0]);
419+ a30 = (e0 - e[0]) / (e[3] - e[0]);
420+
421+ *v = 3.0 * a10 * a20 / (e[3] - e[0]);
422+
423+ tsmall[0][0] = 1.0 - a10;
424+ tsmall[0][1] = 1.0 - a20;
425+ tsmall[0][2] = 1.0 - a30;
426+ tsmall[1][0] = a10;
427+ tsmall[1][1] = 0.0;
428+ tsmall[1][2] = 0.0;
429+ tsmall[2][0] = 0.0;
430+ tsmall[2][1] = a20;
431+ tsmall[2][2] = 0.0;
432+ tsmall[3][0] = 0.0;
433+ tsmall[3][1] = 0.0;
434+ tsmall[3][2] = a30;
435+}
436+/*
437+Cut triangle B1
438+*/
439+void libtetrabz_triangle_b1(
440+ double *e,
441+ double e0,
442+ double *v,
443+ double tsmall[4][3]
444+) {
445+ double a30, a13, a20;
446+ a30 = (e0 - e[0]) / (e[3] - e[0]);
447+ a13 = (e0 - e[3]) / (e[1] - e[3]);
448+ a20 = (e0 - e[0]) / (e[2] - e[0]);
449+
450+ *v = 3.0 * a30 * a13 / (e[2] - e[0]);
451+
452+ tsmall[0][0] = 1.0 - a20;
453+ tsmall[0][1] = 1.0 - a30;
454+ tsmall[0][2] = 0.0;
455+ tsmall[1][0] = 0.0;
456+ tsmall[1][1] = 0.0;
457+ tsmall[1][2] = a13;
458+ tsmall[2][0] = a20;
459+ tsmall[2][1] = 0.0;
460+ tsmall[2][2] = 0.0;
461+ tsmall[3][0] = 0.0;
462+ tsmall[3][1] = a30;
463+ tsmall[3][2] = 1.0 - a13;
464+}
465+/*
466+Cut triangle B2
467+*/
468+void libtetrabz_triangle_b2(
469+ double *e,
470+ double e0,
471+ double *v,
472+ double tsmall[4][3]
473+)
474+{
475+ double a12, a31, a20;
476+ a12 = (e0 - e[2]) / (e[1] - e[2]);
477+ a31 = (e0 - e[1]) / (e[3] - e[1]);
478+ a20 = (e0 - e[0]) / (e[2] - e[0]);
479+
480+ *v = 3.0 * a12 * a31 / (e[2] - e[0]);
481+
482+ tsmall[0][0] = 1.0 - a20;
483+ tsmall[0][1] = 0.0;
484+ tsmall[0][2] = 0.0;
485+ tsmall[1][0] = 0.0;
486+ tsmall[1][1] = a12;
487+ tsmall[1][2] = 1.0 - a31;
488+ tsmall[2][0] = a20;
489+ tsmall[2][1] = 1.0 - a12;
490+ tsmall[2][2] = 0.0;
491+ tsmall[3][0] = 0.0;
492+ tsmall[3][1] = 0.0;
493+ tsmall[3][2]= a31;
494+}
495+/*
496+Cut triangle C1
497+*/
498+void libtetrabz_triangle_c1(
499+ double *e,
500+ double e0,
501+ double *v,
502+ double tsmall[4][3]
503+)
504+{
505+ double a03, a13, a23;
506+ a03 = (e0 - e[3]) / (e[0] - e[3]);
507+ a13 = (e0 - e[3]) / (e[1] - e[3]);
508+ a23 = (e0 - e[3]) / (e[2] - e[3]);
509+
510+ *v = 3.0 * a03 * a13 / (e[3] - e[2]);
511+
512+ tsmall[0][0] = a03;
513+ tsmall[0][1] = 0.0;
514+ tsmall[0][2] = 0.0;
515+ tsmall[1][0] = 0.0;
516+ tsmall[1][1] = a13;
517+ tsmall[1][2] = 0.0;
518+ tsmall[2][0] = 0.0;
519+ tsmall[2][1] = 0.0;
520+ tsmall[2][2] = a23;
521+ tsmall[3][0] = 1.0 - a03;
522+ tsmall[3][1] = 1.0 - a13;
523+ tsmall[3][2] = 1.0 - a23;
524+}
525+/*
526+Sort eigenvalues
527+*/
528+void eig_sort(
529+ int n, //!< [in] the number of components
530+ double *key, //!< [in] Variables to be sorted [n].
531+ int *swap //!< [out] Order of index (sorted)
532+)
533+{
534+ int i, j, k, min_loc;
535+ double min_val;
536+
537+ for (i = 0; i < n; ++i) swap[i] = i;
538+
539+ for (i = 0; i < n - 1; ++i) {
540+ min_val = key[i];
541+ min_loc = i;
542+ for (j = i + 1; j < n; ++j) {
543+ if (min_val > key[j]) {
544+ min_val = key[j];
545+ min_loc = j;
546+ }
547+ }
548+ if (key[i] > min_val) {
549+ /*
550+ Swap
551+ */
552+ key[min_loc] = key[i];
553+ key[i] = min_val;
554+
555+ k = swap[min_loc];
556+ swap[min_loc] = swap[i];
557+ swap[i] = k;
558+ }
559+ }/*for (i = 0; i < n - 1; ++i)*/
560+}/*eig_sort*/
561+/*
562+2nd step of tetrahedron method.
563+*/
564+static void libtetrabz_dbldelta2(
565+ int nb,
566+ double **ej,
567+ double **w
568+) {
569+ int i3, ib, indx[3];
570+ double a10, a20, a02, a12, v, e[3], e_abs;
571+
572+ for (ib = 0; ib < nb; ib++) {
573+
574+ e_abs = 0.0;
575+ for (i3 = 0; i3 < 3; i3++) {
576+ e[i3] = ej[ib][i3];
577+ if (e_abs < fabs(e[i3])) e_abs = fabs(e[i3]);
578+ }
579+ eig_sort(3, e, indx);
580+
581+ if (e_abs < 1.0e-10) {
582+ printf("Nesting ##\n");
583+ }
584+
585+ if ((e[0] < 0.0 && 0.0 <= e[1]) || (e[0] <= 0.0 && 0.0 < e[1])) {
586+
587+ a10 = (0.0 - e[0]) / (e[1] - e[0]);
588+ a20 = (0.0 - e[0]) / (e[2] - e[0]);
589+
590+ v = a10 / (e[2] - e[0]);
591+
592+ w[ib][indx[0]] = v * (2.0 - a10 - a20);
593+ w[ib][indx[1]] = v * a10;
594+ w[ib][indx[2]] = v * a20;
595+ }
596+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
597+
598+ a02 = (0.0 - e[2]) / (e[0] - e[2]);
599+ a12 = (0.0 - e[2]) / (e[1] - e[2]);
600+
601+ v = a12 / (e[2] - e[0]);
602+
603+ w[ib][indx[0]] = v * a02;
604+ w[ib][indx[1]] = v * a12;
605+ w[ib][indx[2]] = v * (2.0 - a02 - a12);
606+ }
607+ else {
608+ for (i3 = 0; i3 < 3; i3++)
609+ w[ib][i3] = 0.0;
610+ }
611+ }
612+}
613+/*
614+Main SUBROUTINE for Delta(E1) * Delta(E2)
615+*/
616+static PyObject* dbldelta(PyObject* self, PyObject* args)
617+{
618+ int it, ik, ib, i20, i4, j3, jb, ** ikv, indx[4], ierr, ng[3], nk, nb;
619+ double wlsm[20][4], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[4][3], thr,
620+ bvec[3][3], ** eig1, ** eig2, *** wght;
621+ PyObject* eig1_po, * eig2_po, * wght_po;
622+ /*
623+ Read input from python object
624+ */
625+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
626+ &ng[0], &ng[1], &ng[2], &nk, &nb,
627+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
628+ &eig1_po, &eig2_po))
629+ return NULL;
630+ /*
631+ convert python list object to array
632+ */
633+ eig1 = (double**)malloc(nk * sizeof(double*));
634+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
635+ eig2 = (double**)malloc(nk * sizeof(double*));
636+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
637+ wght = (double***)malloc(nk * sizeof(double**));
638+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
639+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
640+ for (ik = 0; ik < nk; ik++) {
641+ eig1[ik] = eig1[0] + ik * nb;
642+ eig2[ik] = eig2[0] + ik * nb;
643+ wght[ik] = wght[0] + ik * nb;
644+ for (ib = 0; ib < nb; ib++) {
645+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
646+ }
647+ }
648+
649+ for (ik = 0; ik < nk; ik++)
650+ for (ib = 0; ib < nb; ib++) {
651+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
652+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
653+ }
654+ /*
655+ Start main calculation
656+ */
657+ thr = 1.0e-10;
658+
659+ ikv = (int**)malloc(6 * nk * sizeof(int*));
660+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
661+ for (ik = 0; ik < 6 * nk; ik++) {
662+ ikv[ik] = ikv[0] + ik * 20;
663+ }
664+
665+ ei1 = (double**)malloc(4 * sizeof(double*));
666+ ej1 = (double**)malloc(4 * sizeof(double*));
667+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
668+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
669+ for (i4 = 0; i4 < 4; i4++) {
670+ ei1[i4] = ei1[0] + i4 * nb;
671+ ej1[i4] = ej1[0] + i4 * nb;
672+ }
673+
674+ w1 = (double***)malloc(nb * sizeof(double**));
675+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
676+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
677+ for (ib = 0; ib < nb; ib++) {
678+ w1[ib] = w1[0] + ib * 4;
679+ for (i4 = 0; i4 < 4; i4++) {
680+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
681+ }
682+ }
683+
684+ ej2 = (double**)malloc(nb * sizeof(double*));
685+ ej2[0] = (double*)malloc(nb * 3 * sizeof(double));
686+ for (ib = 0; ib < nb; ib++) {
687+ ej2[ib] = ej2[0] + ib * 3;
688+ }
689+
690+ w2 = (double**)malloc(nb * sizeof(double*));
691+ w2[0] = (double*)malloc(nb * 3 * sizeof(double));
692+ for (ib = 0; ib < nb; ib++) {
693+ w2[ib] = w2[0] + ib * 3;
694+ }
695+
696+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
697+
698+ for (ik = 0; ik < nk; ik++)
699+ for (ib = 0; ib < nb; ib++)
700+ for (jb = 0; jb < nb; jb++)
701+ wght[ik][ib][jb] = 0.0;
702+
703+ for (it = 0; it < 6 * nk; it++) {
704+
705+ for (i4 = 0; i4 < 4; i4++)
706+ for (ib = 0; ib < nb; ib++) {
707+ ei1[i4][ib] = 0.0;
708+ ej1[i4][ib] = 0.0;
709+ }
710+ for (i20 = 0; i20 < 20; i20++) {
711+ for (i4 = 0; i4 < 4; i4++) {
712+ for (ib = 0; ib < nb; ib++) {
713+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
714+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
715+ }
716+ }
717+ }
718+
719+ for (ib = 0; ib < nb; ib++) {
720+
721+ for (i4 = 0; i4 < 4; i4++)
722+ for (jb = 0; jb < nb; jb++)
723+ w1[ib][i4][jb] = 0.0;
724+
725+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
726+ eig_sort(4, e, indx);
727+
728+ if (e[0] < 0.0 && 0.0 <= e[1]) {
729+
730+ libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
731+
732+ if (v > thr) {
733+ for (jb = 0; jb < nb; jb++)
734+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
735+ for (i4 = 0; i4 < 4; i4++)
736+ for (jb = 0; jb < nb; jb++)
737+ for (j3 = 0; j3 < 3; j3++)
738+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
739+ libtetrabz_dbldelta2(nb, ej2, w2);
740+ for (i4 = 0; i4 < 4; i4++)
741+ for (jb = 0; jb < nb; jb++)
742+ for (j3 = 0; j3 < 3; j3++)
743+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
744+ }
745+ }
746+ else if (e[1] < 0.0 && 0.0 <= e[2]) {
747+
748+ libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
749+
750+ if (v > thr) {
751+ for (jb = 0; jb < nb; jb++)
752+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
753+ for (i4 = 0; i4 < 4; i4++)
754+ for (jb = 0; jb < nb; jb++)
755+ for (j3 = 0; j3 < 3; j3++)
756+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
757+ libtetrabz_dbldelta2(nb, ej2, w2);
758+ for (i4 = 0; i4 < 4; i4++)
759+ for (jb = 0; jb < nb; jb++)
760+ for (j3 = 0; j3 < 3; j3++)
761+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
762+ }
763+
764+ libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
765+
766+ if (v > thr) {
767+ for (jb = 0; jb < nb; jb++)
768+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
769+ for (i4 = 0; i4 < 4; i4++)
770+ for (jb = 0; jb < nb; jb++)
771+ for (j3 = 0; j3 < 3; j3++)
772+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
773+ libtetrabz_dbldelta2(nb, ej2, w2);
774+ for (i4 = 0; i4 < 4; i4++)
775+ for (jb = 0; jb < nb; jb++)
776+ for (j3 = 0; j3 < 3; j3++)
777+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
778+ }
779+ }
780+ else if (e[2] < 0.0 && 0.0 < e[3]) {
781+
782+ libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
783+
784+ if (v > thr) {
785+ for (jb = 0; jb < nb; jb++)
786+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
787+ for (i4 = 0; i4 < 4; i4++)
788+ for (jb = 0; jb < nb; jb++)
789+ for (j3 = 0; j3 < 3; j3++)
790+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
791+ libtetrabz_dbldelta2(nb, ej2, w2);
792+ for (i4 = 0; i4 < 4; i4++)
793+ for (jb = 0; jb < nb; jb++)
794+ for (j3 = 0; j3 < 3; j3++)
795+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
796+ }
797+ }
798+ else {
799+ continue;
800+ }
801+ }
802+ for (i20 = 0; i20 < 20; i20++)
803+ for (ib = 0; ib < nb; ib++)
804+ for (i4 = 0; i4 < 4; i4++)
805+ for (jb = 0; jb < nb; jb++)
806+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
807+ }
808+ for (ik = 0; ik < nk; ik++)
809+ for (ib = 0; ib < nb; ib++)
810+ for (jb = 0; jb < nb; jb++)
811+ wght[ik][ib][jb] /= (6.0 * (double)nk);
812+
813+ free(ikv[0]);
814+ free(ikv);
815+ free(ei1[0]);
816+ free(ei1);
817+ free(ej1[0]);
818+ free(ej1);
819+ free(ej2[0]);
820+ free(ej2);
821+ free(w1[0][0]);
822+ free(w1[0]);
823+ free(w1);
824+ free(w2[0]);
825+ free(w2);
826+ /*
827+ Convert weight to python list object
828+ */
829+ wght_po = PyList_New(nk * nb * nb);
830+ ierr = 0;
831+ for (ik = 0; ik < nk; ik++)
832+ for (ib = 0; ib < nb; ib++)
833+ for (jb = 0; jb < nb; jb++)
834+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
835+ if (ierr != 0)printf("Error in PyList_SetItem\n");
836+
837+ free(eig1[0]);
838+ free(eig1);
839+ free(eig2[0]);
840+ free(eig2);
841+ free(wght[0][0]);
842+ free(wght[0]);
843+ free(wght);
844+
845+ return wght_po;
846+}
847+/*
848+Tetrahedron method for theta( - de)
849+*/
850+static void libtetrabz_dblstep2(
851+ int nb,
852+ double ei[4],
853+ double **ej,
854+ double **w1
855+) {
856+ int i4, j4, ib, indx[4];
857+ double v, thr, e[4], tsmall[4][4];
858+
859+ thr = 1.0e-8;
860+
861+ for (ib = 0; ib < nb; ib++) {
862+
863+ for (i4 = 0; i4 < 4; i4++)
864+ w1[ib][i4] = 0.0;
865+
866+ for (i4 = 0; i4 < 4; i4++) e[i4] = - ei[i4] + ej[ib][i4];
867+ eig_sort(4, e, indx);
868+
869+ if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
870+ /*
871+ Theta(0) = 0.5
872+ */
873+ v = 0.5;
874+ for (i4 = 0; i4 < 4; i4++)
875+ w1[ib][i4] += v * 0.25;
876+ }
877+ else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
878+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
879+ for (i4 = 0; i4 < 4; i4++)
880+ for (j4 = 0; j4 < 4; j4++)
881+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
882+ }
883+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
884+
885+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
886+ for (i4 = 0; i4 < 4; i4++)
887+ for (j4 = 0; j4 < 4; j4++)
888+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
889+
890+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
891+ for (i4 = 0; i4 < 4; i4++)
892+ for (j4 = 0; j4 < 4; j4++)
893+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
894+
895+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
896+ for (i4 = 0; i4 < 4; i4++)
897+ for (j4 = 0; j4 < 4; j4++)
898+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
899+ }
900+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
901+
902+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
903+ for (i4 = 0; i4 < 4; i4++)
904+ for (j4 = 0; j4 < 4; j4++)
905+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
906+
907+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
908+ for (i4 = 0; i4 < 4; i4++)
909+ for (j4 = 0; j4 < 4; j4++)
910+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
911+
912+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
913+ for (i4 = 0; i4 < 4; i4++)
914+ for (j4 = 0; j4 < 4; j4++)
915+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
916+ }
917+ else if (e[3] <= 0.0) {
918+ for (i4 = 0; i4 < 4; i4++)
919+ w1[ib][i4] += 0.25;
920+ }
921+ }
922+}
923+/*
924+Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
925+*/
926+static PyObject* dblstep(PyObject* self, PyObject* args)
927+{
928+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4], ierr, ng[3], nk, nb;
929+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr,
930+ bvec[3][3], ** eig1, ** eig2, *** wght;
931+ PyObject* eig1_po, * eig2_po, * wght_po;
932+
933+ /*
934+ Read input from python object
935+ */
936+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
937+ &ng[0], &ng[1], &ng[2], &nk, &nb,
938+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
939+ &eig1_po, &eig2_po))
940+ return NULL;
941+ /*
942+ convert python list object to array
943+ */
944+ eig1 = (double**)malloc(nk * sizeof(double*));
945+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
946+ eig2 = (double**)malloc(nk * sizeof(double*));
947+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
948+ wght = (double***)malloc(nk * sizeof(double**));
949+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
950+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
951+ for (ik = 0; ik < nk; ik++) {
952+ eig1[ik] = eig1[0] + ik * nb;
953+ eig2[ik] = eig2[0] + ik * nb;
954+ wght[ik] = wght[0] + ik * nb;
955+ for (ib = 0; ib < nb; ib++) {
956+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
957+ }
958+ }
959+
960+ for (ik = 0; ik < nk; ik++)
961+ for (ib = 0; ib < nb; ib++) {
962+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
963+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
964+ }
965+ /*
966+ Start main calculation
967+ */
968+ ikv = (int**)malloc(6 * nk * sizeof(int*));
969+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
970+ for (ik = 0; ik < 6 * nk; ik++) {
971+ ikv[ik] = ikv[0] + ik * 20;
972+ }
973+
974+ ei1 = (double**)malloc(4 * sizeof(double*));
975+ ej1 = (double**)malloc(4 * sizeof(double*));
976+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
977+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
978+ for (i4 = 0; i4 < 4; i4++) {
979+ ei1[i4] = ei1[0] + i4 * nb;
980+ ej1[i4] = ej1[0] + i4 * nb;
981+ }
982+
983+ ej2 = (double**)malloc(nb * sizeof(double*));
984+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
985+ for (ib = 0; ib < nb; ib++) {
986+ ej2[ib] = ej2[0] + ib * 4;
987+ }
988+
989+ w1 = (double***)malloc(nb * sizeof(double**));
990+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
991+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
992+ for (ib = 0; ib < nb; ib++) {
993+ w1[ib] = w1[0] + ib * 4;
994+ for (i4 = 0; i4 < 4; i4++) {
995+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
996+ }
997+ }
998+
999+ w2 = (double**)malloc(nb * sizeof(double*));
1000+ w2[0] = (double*)malloc(nb * 4 * sizeof(double));
1001+ for (ib = 0; ib < nb; ib++) {
1002+ w2[ib] = w2[0] + ib * 4;
1003+ }
1004+
1005+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
1006+
1007+ for (ik = 0; ik < nk; ik++)
1008+ for (ib = 0; ib < nb; ib++)
1009+ for (jb = 0; jb < nb; jb++)
1010+ wght[ik][ib][jb] = 0.0;
1011+
1012+ thr = 1.0e-8;
1013+
1014+ for(it = 0; it < 6*nk; it++) {
1015+
1016+ for (i4 = 0; i4 < 4; i4++)
1017+ for (ib = 0; ib < nb; ib++) {
1018+ ei1[i4][ib] = 0.0;
1019+ ej1[i4][ib] = 0.0;
1020+ }
1021+ for (i20 = 0; i20 < 20; i20++) {
1022+ for (i4 = 0; i4 < 4; i4++) {
1023+ for (ib = 0; ib < nb; ib++) {
1024+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
1025+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
1026+ }
1027+ }
1028+ }
1029+
1030+ for (ib = 0; ib < nb; ib++) {
1031+
1032+ for (i4 = 0; i4 < 4; i4++)
1033+ for (jb = 0; jb < nb; jb++)
1034+ w1[ib][i4][jb] = 0.0;
1035+
1036+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
1037+ eig_sort(4, e, indx);
1038+
1039+ if (e[0] <= 0.0 && 0.0 < e[1]) {
1040+
1041+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
1042+
1043+ if (v > thr) {
1044+ for (j4 = 0; j4 < 4; j4++) {
1045+ ei2[j4] = 0.0;
1046+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1047+ }
1048+ for (i4 = 0; i4 < 4; i4++)
1049+ for (j4 = 0; j4 < 4; j4++) {
1050+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1051+ for (jb = 0; jb < nb; jb++)
1052+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1053+ }
1054+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1055+ for (i4 = 0; i4 < 4; i4++)
1056+ for (jb = 0; jb < nb; jb++)
1057+ for (j4 = 0; j4 < 4; j4++)
1058+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1059+ }
1060+ }
1061+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
1062+
1063+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
1064+
1065+ if (v > thr) {
1066+ for (j4 = 0; j4 < 4; j4++) {
1067+ ei2[j4] = 0.0;
1068+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1069+ }
1070+ for (i4 = 0; i4 < 4; i4++)
1071+ for (j4 = 0; j4 < 4; j4++) {
1072+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1073+ for (jb = 0; jb < nb; jb++)
1074+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1075+ }
1076+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1077+ for (i4 = 0; i4 < 4; i4++)
1078+ for (jb = 0; jb < nb; jb++)
1079+ for (j4 = 0; j4 < 4; j4++)
1080+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1081+ }
1082+
1083+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
1084+
1085+ if (v > thr) {
1086+ for (j4 = 0; j4 < 4; j4++) {
1087+ ei2[j4] = 0.0;
1088+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1089+ }
1090+ for (i4 = 0; i4 < 4; i4++)
1091+ for (j4 = 0; j4 < 4; j4++) {
1092+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1093+ for (jb = 0; jb < nb; jb++)
1094+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1095+ }
1096+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1097+ for (i4 = 0; i4 < 4; i4++)
1098+ for (jb = 0; jb < nb; jb++)
1099+ for (j4 = 0; j4 < 4; j4++)
1100+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1101+ }
1102+
1103+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
1104+
1105+ if (v > thr) {
1106+ for (j4 = 0; j4 < 4; j4++) {
1107+ ei2[j4] = 0.0;
1108+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1109+ }
1110+ for (i4 = 0; i4 < 4; i4++)
1111+ for (j4 = 0; j4 < 4; j4++) {
1112+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1113+ for (jb = 0; jb < nb; jb++)
1114+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1115+ }
1116+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1117+ for (i4 = 0; i4 < 4; i4++)
1118+ for (jb = 0; jb < nb; jb++)
1119+ for (j4 = 0; j4 < 4; j4++)
1120+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1121+ }
1122+ }
1123+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
1124+
1125+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
1126+
1127+ if (v > thr) {
1128+ for (j4 = 0; j4 < 4; j4++) {
1129+ ei2[j4] = 0.0;
1130+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1131+ }
1132+ for (i4 = 0; i4 < 4; i4++)
1133+ for (j4 = 0; j4 < 4; j4++) {
1134+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1135+ for (jb = 0; jb < nb; jb++)
1136+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1137+ }
1138+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1139+ for (i4 = 0; i4 < 4; i4++)
1140+ for (jb = 0; jb < nb; jb++)
1141+ for (j4 = 0; j4 < 4; j4++)
1142+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1143+ }
1144+
1145+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
1146+
1147+ if (v > thr) {
1148+ for (j4 = 0; j4 < 4; j4++) {
1149+ ei2[j4] = 0.0;
1150+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1151+ }
1152+ for (i4 = 0; i4 < 4; i4++)
1153+ for (j4 = 0; j4 < 4; j4++) {
1154+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1155+ for (jb = 0; jb < nb; jb++)
1156+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1157+ }
1158+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1159+ for (i4 = 0; i4 < 4; i4++)
1160+ for (jb = 0; jb < nb; jb++)
1161+ for (j4 = 0; j4 < 4; j4++)
1162+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1163+ }
1164+
1165+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
1166+
1167+ if (v > thr) {
1168+ for (j4 = 0; j4 < 4; j4++) {
1169+ ei2[j4] = 0.0;
1170+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1171+ }
1172+ for (i4 = 0; i4 < 4; i4++)
1173+ for (j4 = 0; j4 < 4; j4++) {
1174+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1175+ for (jb = 0; jb < nb; jb++)
1176+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1177+ }
1178+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1179+ for (i4 = 0; i4 < 4; i4++)
1180+ for (jb = 0; jb < nb; jb++)
1181+ for (j4 = 0; j4 < 4; j4++)
1182+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
1183+ }
1184+ }
1185+ else if (e[3] <= 0.0) {
1186+ for (i4 = 0; i4 < 4; i4++) {
1187+ ei2[i4] = ei1[i4][ib];
1188+ for (jb = 0; jb < nb; jb++)
1189+ ej2[jb][i4] = ej1[i4][jb];
1190+ }
1191+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
1192+ for (i4 = 0; i4 < 4; i4++)
1193+ for (jb = 0; jb < nb; jb++)
1194+ w1[ib][i4][jb] += w2[jb][i4];
1195+ }
1196+ else {
1197+ continue;
1198+ }
1199+ }
1200+ for (i20 = 0; i20 < 20; i20++)
1201+ for (ib = 0; ib < nb; ib++)
1202+ for (i4 = 0; i4 < 4; i4++)
1203+ for (jb = 0; jb < nb; jb++)
1204+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
1205+ }
1206+ for (ik = 0; ik < nk; ik++)
1207+ for (ib = 0; ib < nb; ib++)
1208+ for (jb = 0; jb < nb; jb++)
1209+ wght[ik][ib][jb] /= (6.0 * (double) nk);
1210+
1211+ free(ikv[0]);
1212+ free(ikv);
1213+ free(ei1[0]);
1214+ free(ei1);
1215+ free(ej1[0]);
1216+ free(ej1);
1217+ free(ej2[0]);
1218+ free(ej2);
1219+ free(w1[0][0]);
1220+ free(w1[0]);
1221+ free(w1);
1222+ free(w2[0]);
1223+ free(w2);
1224+ /*
1225+ Convert weight to python list object
1226+ */
1227+ wght_po = PyList_New(nk * nb * nb);
1228+ ierr = 0;
1229+ for (ik = 0; ik < nk; ik++)
1230+ for (ib = 0; ib < nb; ib++)
1231+ for (jb = 0; jb < nb; jb++)
1232+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
1233+ if (ierr != 0)printf("Error in PyList_SetItem\n");
1234+
1235+ free(eig1[0]);
1236+ free(eig1);
1237+ free(eig2[0]);
1238+ free(eig2);
1239+ free(wght[0][0]);
1240+ free(wght[0]);
1241+ free(wght);
1242+
1243+ return wght_po;
1244+}
1245+/*
1246+Compute Dos : Delta(E - E1)
1247+*/
1248+static PyObject* dos(PyObject* self, PyObject* args) {
1249+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4], ierr, ng[3], nk, nb, ne;
1250+ double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][3], bvec[3][3], ** eig, * e0, *** wght;
1251+ PyObject* eig_po, * e0_po, * wght_po;
1252+ /*
1253+ Read input from python object
1254+ */
1255+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOO",
1256+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
1257+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
1258+ &eig_po, &e0_po))
1259+ return NULL;
1260+ /*
1261+ convert python list object to array
1262+ */
1263+ eig = (double**)malloc(nk * sizeof(double*));
1264+ eig[0] = (double*)malloc(nk * nb * sizeof(double));
1265+ wght = (double***)malloc(nk * sizeof(double**));
1266+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
1267+ wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));
1268+ for (ik = 0; ik < nk; ik++) {
1269+ eig[ik] = eig[0] + ik * nb;
1270+ wght[ik] = wght[0] + ik * nb;
1271+ for (ib = 0; ib < nb; ib++) {
1272+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
1273+ }
1274+ }
1275+ e0 = (double*)malloc(ne * sizeof(double));
1276+
1277+ for (ik = 0; ik < nk; ik++)
1278+ for (ib = 0; ib < nb; ib++)
1279+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
1280+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
1281+ /*
1282+ Start main calculation
1283+ */
1284+ ikv = (int**)malloc(6 * nk * sizeof(int*));
1285+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
1286+ for (ik = 0; ik < 6 * nk; ik++) {
1287+ ikv[ik] = ikv[0] + ik * 20;
1288+ }
1289+
1290+ ei1 = (double**)malloc(4 * sizeof(double*));
1291+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
1292+ for (ii = 0; ii < 4; ii++) {
1293+ ei1[ii] = ei1[0] + ii * nb;
1294+ }
1295+
1296+ w1 = (double***)malloc(nb * sizeof(double**));
1297+ w1[0] = (double**)malloc(nb * ne * sizeof(double*));
1298+ w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
1299+ for (ib = 0; ib < nb; ib++) {
1300+ w1[ib] = w1[0] + ib * ne;
1301+ for (ie = 0; ie < ne; ie++) {
1302+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
1303+ }
1304+ }
1305+
1306+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
1307+
1308+ for (ik = 0; ik < nk; ik++)
1309+ for (ib = 0; ib < nb; ib++)
1310+ for (ie = 0; ie < ne; ie++)
1311+ wght[ik][ib][ie] = 0.0;
1312+
1313+ for (it = 0; it < 6 * nk; it++) {
1314+
1315+ for (ii = 0; ii < 4; ii++)
1316+ for (ib = 0; ib < nb; ib++)
1317+ ei1[ii][ib] = 0.0;
1318+ for (jj = 0; jj < 20; jj++)
1319+ for (ii = 0; ii < 4; ii++)
1320+ for (ib = 0; ib < nb; ib++)
1321+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
1322+
1323+ for (ib = 0; ib < nb; ib++) {
1324+
1325+ for (ie = 0; ie < ne; ie++)
1326+ for (ii = 0; ii < 4; ii++)
1327+ w1[ib][ie][ii] = 0.0;
1328+
1329+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
1330+ eig_sort(4, e, indx);
1331+
1332+ for (ie = 0; ie < ne; ie++) {
1333+
1334+ if ((e[0] < e0[ie] && e0[ie] <= e[1]) || (e[0] <= e0[ie] && e0[ie] < e[1])) {
1335+
1336+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
1337+ for (ii = 0; ii < 4; ii++)
1338+ for (jj = 0; jj < 3; jj++)
1339+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
1340+
1341+ }
1342+ else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
1343+
1344+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
1345+ for (ii = 0; ii < 4; ii++)
1346+ for (jj = 0; jj < 3; jj++)
1347+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
1348+
1349+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
1350+ for (ii = 0; ii < 4; ii++)
1351+ for (jj = 0; jj < 3; jj++)
1352+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
1353+ }
1354+ else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
1355+
1356+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
1357+ for (ii = 0; ii < 4; ii++)
1358+ for (jj = 0; jj < 3; jj++)
1359+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
1360+ }
1361+ else {
1362+ continue;
1363+ }
1364+ }
1365+ }
1366+ for (ii = 0; ii < 20; ii++)
1367+ for (ib = 0; ib < nb; ib++)
1368+ for (ie = 0; ie < ne; ie++)
1369+ for (jj = 0; jj < 4; jj++)
1370+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
1371+ }
1372+ for (ik = 0; ik < nk; ik++)
1373+ for (ib = 0; ib < nb; ib++)
1374+ for (ie = 0; ie < ne; ie++)
1375+ wght[ik][ib][ie] /= (6.0 * (double)nk);
1376+
1377+ free(ikv[0]);
1378+ free(ikv);
1379+ free(ei1[0]);
1380+ free(ei1);
1381+ free(w1[0][0]);
1382+ free(w1[0]);
1383+ free(w1);
1384+ /*
1385+ Convert weight to python list object
1386+ */
1387+ wght_po = PyList_New(nk * nb * ne);
1388+ ierr = 0;
1389+ for (ik = 0; ik < nk; ik++)
1390+ for (ib = 0; ib < nb; ib++)
1391+ for (ie = 0; ie < ne; ie++)
1392+ ierr = PyList_SetItem(wght_po, ik * nb * ne + ib * ne + ie, PyFloat_FromDouble(wght[ik][ib][ie]));
1393+ if (ierr != 0)printf("Error in PyList_SetItem\n");
1394+
1395+ free(eig[0]);
1396+ free(eig);
1397+ free(wght[0][0]);
1398+ free(wght[0]);
1399+ free(wght);
1400+ free(e0);
1401+
1402+ return wght_po;
1403+}
1404+/*
1405+Compute integrated Dos : theta(E - E1)
1406+*/
1407+static PyObject* intdos(PyObject* self, PyObject* args) {
1408+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4], ierr, ng[3], nk, nb, ne;
1409+ double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][4], bvec[3][3], ** eig, * e0, *** wght;
1410+ PyObject* eig_po, * e0_po, * wght_po;
1411+ /*
1412+ Read input from python object
1413+ */
1414+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOO",
1415+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
1416+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
1417+ &eig_po, &e0_po))
1418+ return NULL;
1419+ /*
1420+ convert python list object to array
1421+ */
1422+ eig = (double**)malloc(nk * sizeof(double*));
1423+ eig[0] = (double*)malloc(nk * nb * sizeof(double));
1424+ wght = (double***)malloc(nk * sizeof(double**));
1425+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
1426+ wght[0][0] = (double*)malloc(nk * nb * ne * sizeof(double));
1427+ for (ik = 0; ik < nk; ik++) {
1428+ eig[ik] = eig[0] + ik * nb;
1429+ wght[ik] = wght[0] + ik * nb;
1430+ for (ib = 0; ib < nb; ib++) {
1431+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
1432+ }
1433+ }
1434+ e0 = (double*)malloc(ne * sizeof(double));
1435+
1436+ for (ik = 0; ik < nk; ik++)
1437+ for (ib = 0; ib < nb; ib++)
1438+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
1439+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
1440+ /*
1441+ Start main calculation
1442+ */
1443+ ikv = (int**)malloc(6 * nk * sizeof(int*));
1444+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
1445+ for (ik = 0; ik < 6 * nk; ik++) {
1446+ ikv[ik] = ikv[0] + ik * 20;
1447+ }
1448+
1449+ ei1 = (double**)malloc(4 * sizeof(double*));
1450+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
1451+ for (ii = 0; ii < 4; ii++) {
1452+ ei1[ii] = ei1[0] + ii * nb;
1453+ }
1454+
1455+ w1 = (double***)malloc(nb * sizeof(double**));
1456+ w1[0] = (double**)malloc(nb * ne * sizeof(double*));
1457+ w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
1458+ for (ib = 0; ib < nb; ib++) {
1459+ w1[ib] = w1[0] + ib * ne;
1460+ for (ie = 0; ie < ne; ie++) {
1461+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
1462+ }
1463+ }
1464+
1465+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
1466+
1467+ for (ik = 0; ik < nk; ik++)
1468+ for (ib = 0; ib < nb; ib++)
1469+ for (ie = 0; ie < ne; ie++)
1470+ wght[ik][ib][ie] = 0.0;
1471+
1472+ for (it = 0; it < 6 * nk; it++) {
1473+
1474+ for (ii = 0; ii < 4; ii++)
1475+ for (ib = 0; ib < nb; ib++)
1476+ ei1[ii][ib] = 0.0;
1477+ for (jj = 0; jj < 20; jj++)
1478+ for (ii = 0; ii < 4; ii++)
1479+ for (ib = 0; ib < nb; ib++)
1480+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
1481+
1482+ for (ib = 0; ib < nb; ib++)
1483+ for (ie = 0; ie < ne; ie++)
1484+ for (ii = 0; ii < 4; ii++)
1485+ w1[ib][ie][ii] = 0.0;
1486+
1487+ for (ib = 0; ib < nb; ib++) {
1488+
1489+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
1490+ eig_sort(4, e, indx);
1491+
1492+ for (ie = 0; ie < ne; ie++) {
1493+
1494+ if ((e[0] <= e0[ie] && e0[ie] < e[1]) || (e[0] < e0[ie] && e0[ie] <= e[1])) {
1495+ libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
1496+ for (ii = 0; ii < 4; ii++)
1497+ for (jj = 0; jj < 4; jj++)
1498+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1499+ }
1500+ else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
1501+
1502+ libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
1503+ for (ii = 0; ii < 4; ii++)
1504+ for (jj = 0; jj < 4; jj++)
1505+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1506+
1507+ libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
1508+ for (ii = 0; ii < 4; ii++)
1509+ for (jj = 0; jj < 4; jj++)
1510+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1511+
1512+ libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
1513+ for (ii = 0; ii < 4; ii++)
1514+ for (jj = 0; jj < 4; jj++)
1515+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1516+
1517+ }
1518+ else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
1519+
1520+ libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
1521+ for (ii = 0; ii < 4; ii++)
1522+ for (jj = 0; jj < 4; jj++)
1523+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1524+
1525+ libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
1526+ for (ii = 0; ii < 4; ii++)
1527+ for (jj = 0; jj < 4; jj++)
1528+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1529+
1530+ libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
1531+ for (ii = 0; ii < 4; ii++)
1532+ for (jj = 0; jj < 4; jj++)
1533+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
1534+
1535+ }
1536+ else if (e[3] <= e0[ie]) {
1537+ for (ii = 0; ii < 4; ii++)
1538+ w1[ib][ie][ii] = 0.25;
1539+ }
1540+ else {
1541+ continue;
1542+ }
1543+ }
1544+ }
1545+
1546+ for (ii = 0; ii < 20; ii++)
1547+ for (ib = 0; ib < nb; ib++)
1548+ for (ie = 0; ie < ne; ie++)
1549+ for (jj = 0; jj < 4; jj++)
1550+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
1551+ }
1552+ for (ik = 0; ik < nk; ik++)
1553+ for (ib = 0; ib < nb; ib++)
1554+ for (ie = 0; ie < ne; ie++)
1555+ wght[ik][ib][ie] /= (6.0 * (double)nk);
1556+
1557+ free(ikv[0]);
1558+ free(ikv);
1559+ free(ei1[0]);
1560+ free(ei1);
1561+ free(w1[0][0]);
1562+ free(w1[0]);
1563+ free(w1);
1564+ /*
1565+ Convert weight to python list object
1566+ */
1567+ wght_po = PyList_New(nk * nb * ne);
1568+ ierr = 0;
1569+ for (ik = 0; ik < nk; ik++)
1570+ for (ib = 0; ib < nb; ib++)
1571+ for (ie = 0; ie < ne; ie++)
1572+ ierr = PyList_SetItem(wght_po, ik * nb * ne + ib * ne + ie, PyFloat_FromDouble(wght[ik][ib][ie]));
1573+ if (ierr != 0)printf("Error in PyList_SetItem\n");
1574+
1575+ free(eig[0]);
1576+ free(eig);
1577+ free(wght[0][0]);
1578+ free(wght[0]);
1579+ free(wght);
1580+ free(e0);
1581+
1582+ return wght_po; }
1583+/*
1584+ 3rd step for Fermi's golden rule
1585+*/
1586+static void libtetrabz_fermigr3(
1587+ int ne,
1588+ double *e0,
1589+ double de[4],
1590+ double **w1
1591+) {
1592+ int i4, j3, ie, indx[4];
1593+ double e[4], tsmall[4][3], v;
1594+
1595+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
1596+ eig_sort(4, e, indx);
1597+
1598+ for (ie = 0; ie < ne; ie++) {
1599+
1600+ for (i4 = 0; i4 < 4; i4++) w1[i4][ie] = 0.0;
1601+
1602+ if (e[0] < e0[ie] && e0[ie] <= e[1]) {
1603+
1604+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
1605+ for (i4 = 0; i4 < 4; i4++)
1606+ for (j3 = 0; j3 < 3; j3++)
1607+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
1608+ }
1609+ else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
1610+
1611+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
1612+ for (i4 = 0; i4 < 4; i4++)
1613+ for (j3 = 0; j3 < 3; j3++)
1614+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
1615+
1616+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
1617+ for (i4 = 0; i4 < 4; i4++)
1618+ for (j3 = 0; j3 < 3; j3++)
1619+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
1620+ }
1621+ else if (e[2] < e0[ie] && e0[ie] < e[3]) {
1622+
1623+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
1624+ for (i4 = 0; i4 < 4; i4++)
1625+ for (j3 = 0; j3 < 3; j3++)
1626+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
1627+ }
1628+ }
1629+}
1630+/*
1631+2nd step for Fermi's golden rule
1632+*/
1633+static void libtetrabz_fermigr2(
1634+ int nb,
1635+ int ne,
1636+ double *e0,
1637+ double *ei1,
1638+ double **ej1,
1639+ double ***w1
1640+) {
1641+ int ib, i4, j4, ie, indx[4];
1642+ double e[4], tsmall[4][4], v, de[4], thr, **w2;
1643+
1644+ w2 = (double**)malloc(4 * sizeof(double*));
1645+ w2[0] = (double*)malloc(4 * ne * sizeof(double));
1646+ for (i4 = 0; i4 < 4; i4++) {
1647+ w2[i4] = w2[0] + i4 * ne;
1648+ }
1649+
1650+ thr = 1.0e-8;
1651+
1652+ for (ib = 0; ib < nb; ib++) {
1653+
1654+ for (i4 = 0; i4 < 4; i4++)
1655+ for (ie = 0; ie < ne; ie++)
1656+ w1[ib][i4][ie] = 0.0;
1657+
1658+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
1659+ eig_sort(4, e, indx);
1660+
1661+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
1662+
1663+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
1664+
1665+ if (v > thr) {
1666+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1667+ for (i4 = 0; i4 < 4; i4++)
1668+ for (j4 = 0; j4 < 4; j4++)
1669+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1670+ libtetrabz_fermigr3(ne, e0, de, w2);
1671+ for (i4 = 0; i4 < 4; i4++)
1672+ for (j4 = 0; j4 < 4; j4++)
1673+ for (ie = 0; ie < ne; ie++)
1674+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1675+ }
1676+ }
1677+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
1678+
1679+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
1680+
1681+ if (v > thr) {
1682+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1683+ for (i4 = 0; i4 < 4; i4++)
1684+ for (j4 = 0; j4 < 4; j4++)
1685+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1686+ libtetrabz_fermigr3(ne, e0, de, w2);
1687+ for (i4 = 0; i4 < 4; i4++)
1688+ for (j4 = 0; j4 < 4; j4++)
1689+ for (ie = 0; ie < ne; ie++)
1690+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1691+ }
1692+
1693+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
1694+
1695+ if (v > thr) {
1696+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1697+ for (i4 = 0; i4 < 4; i4++)
1698+ for (j4 = 0; j4 < 4; j4++)
1699+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1700+ libtetrabz_fermigr3(ne, e0, de, w2);
1701+ for (i4 = 0; i4 < 4; i4++)
1702+ for (j4 = 0; j4 < 4; j4++)
1703+ for (ie = 0; ie < ne; ie++)
1704+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1705+ }
1706+
1707+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
1708+
1709+ if (v > thr) {
1710+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1711+ for (i4 = 0; i4 < 4; i4++)
1712+ for (j4 = 0; j4 < 4; j4++)
1713+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1714+ libtetrabz_fermigr3(ne, e0, de, w2);
1715+ for (i4 = 0; i4 < 4; i4++)
1716+ for (j4 = 0; j4 < 4; j4++)
1717+ for (ie = 0; ie < ne; ie++)
1718+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1719+ }
1720+ }
1721+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
1722+
1723+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
1724+
1725+ if (v > thr) {
1726+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1727+ for (i4 = 0; i4 < 4; i4++)
1728+ for (j4 = 0; j4 < 4; j4++)
1729+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1730+ libtetrabz_fermigr3(ne, e0, de, w2);
1731+ for (i4 = 0; i4 < 4; i4++)
1732+ for (j4 = 0; j4 < 4; j4++)
1733+ for (ie = 0; ie < ne; ie++)
1734+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1735+ }
1736+
1737+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
1738+
1739+ if (v > thr) {
1740+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1741+ for (i4 = 0; i4 < 4; i4++)
1742+ for (j4 = 0; j4 < 4; j4++)
1743+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1744+ libtetrabz_fermigr3(ne, e0, de, w2);
1745+ for (i4 = 0; i4 < 4; i4++)
1746+ for (j4 = 0; j4 < 4; j4++)
1747+ for (ie = 0; ie < ne; ie++)
1748+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1749+ }
1750+
1751+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
1752+
1753+ if (v > thr) {
1754+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
1755+ for (i4 = 0; i4 < 4; i4++)
1756+ for (j4 = 0; j4 < 4; j4++)
1757+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
1758+ libtetrabz_fermigr3(ne, e0, de, w2);
1759+ for (i4 = 0; i4 < 4; i4++)
1760+ for (j4 = 0; j4 < 4; j4++)
1761+ for (ie = 0; ie < ne; ie++)
1762+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
1763+ }
1764+ }
1765+ else if (e[3] <= 0.0) {
1766+ for (i4 = 0; i4 < 4; i4++)
1767+ de[i4] = ej1[ib][i4] - ei1[i4];
1768+ libtetrabz_fermigr3(ne, e0, de, w2);
1769+ for (i4 = 0; i4 < 4; i4++)
1770+ for (ie = 0; ie < ne; ie++)
1771+ w1[ib][i4][ie] += w2[i4][ie];
1772+ }
1773+ }
1774+
1775+ free(w2[0]);
1776+ free(w2);
1777+}
1778+/*
1779+Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
1780+*/
1781+static PyObject* fermigr(PyObject* self, PyObject* args)
1782+{
1783+ int it, ik, ib, i20, i4, j4, jb, ie, **ikv, indx[4], ierr, ng[3], nk, nb, ne;
1784+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], **** w1, *** w2, v, tsmall[4][4], thr,
1785+ bvec[3][3], ** eig1, ** eig2, *e0, * ***wght;
1786+ PyObject* eig1_po, * eig2_po, *e0_po, * wght_po;
1787+ /*
1788+ Read input from python object
1789+ */
1790+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOOO",
1791+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
1792+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
1793+ &eig1_po, &eig2_po, &e0_po))
1794+ return NULL;
1795+ /*
1796+ convert python list object to array
1797+ */
1798+ eig1 = (double**)malloc(nk * sizeof(double*));
1799+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
1800+ eig2 = (double**)malloc(nk * sizeof(double*));
1801+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
1802+ wght = (double****)malloc(nk * sizeof(double***));
1803+ wght[0] = (double***)malloc(nk * nb * sizeof(double**));
1804+ wght[0][0] = (double**)malloc(nk * nb * nb * sizeof(double*));
1805+ wght[0][0][0] = (double*)malloc(nk * nb * nb *ne *sizeof(double));
1806+ for (ik = 0; ik < nk; ik++) {
1807+ eig1[ik] = eig1[0] + ik * nb;
1808+ eig2[ik] = eig2[0] + ik * nb;
1809+ wght[ik] = wght[0] + ik * nb;
1810+ for (ib = 0; ib < nb; ib++) {
1811+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
1812+ for (jb = 0; jb < nb; jb++) {
1813+ wght[ik][ib][jb] = wght[0][0][0]
1814+ + ik * nb * nb * ne + ib * nb * ne + jb * ne;
1815+ }
1816+ }
1817+ }
1818+ e0 = (double*)malloc(ne * sizeof(double));
1819+
1820+ for (ik = 0; ik < nk; ik++)
1821+ for (ib = 0; ib < nb; ib++) {
1822+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
1823+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
1824+ }
1825+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
1826+ /*
1827+ Start main calculation
1828+ */
1829+ ikv = (int**)malloc(6 * nk * sizeof(int*));
1830+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
1831+ for (it = 0; it < 6 * nk; it++) {
1832+ ikv[it] = ikv[0] + it * 20;
1833+ }
1834+
1835+ ei1 = (double**)malloc(4 * sizeof(double*));
1836+ ej1 = (double**)malloc(4 * sizeof(double*));
1837+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
1838+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
1839+ for (i4 = 0; i4 < 4; i4++) {
1840+ ei1[i4] = ei1[0] + i4 * nb;
1841+ ej1[i4] = ej1[0] + i4 * nb;
1842+ }
1843+
1844+ ej2 = (double**)malloc(nb * sizeof(double*));
1845+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
1846+ for (ib = 0; ib < nb; ib++)
1847+ ej2[ib] = ej2[0] + ib * 4;
1848+
1849+ w1 = (double****)malloc(nb * sizeof(double***));
1850+ w1[0] = (double***)malloc(nb * 4 * sizeof(double**));
1851+ w1[0][0] = (double**)malloc(nb * 4 * nb * sizeof(double*));
1852+ w1[0][0][0] = (double*)malloc(nb * 4 * nb * ne * sizeof(double));
1853+ for (ib = 0; ib < nb; ib++) {
1854+ w1[ib] = w1[0] + ib * 4;
1855+ for (i4 = 0; i4 < 4; i4++) {
1856+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
1857+ for (jb = 0; jb < nb; jb++) {
1858+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
1859+ }
1860+ }
1861+ }
1862+
1863+ w2 = (double***)malloc(nb * sizeof(double**));
1864+ w2[0] = (double**)malloc(nb * 4 * sizeof(double*));
1865+ w2[0][0] = (double*)malloc(nb * 4 * ne * sizeof(double));
1866+ for (ib = 0; ib < nb; ib++) {
1867+ w2[ib] = w2[0] + ib * 4;
1868+ for (i4 = 0; i4 < 4; i4++) {
1869+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
1870+ }
1871+ }
1872+
1873+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
1874+
1875+ for (ik = 0; ik < nk; ik++)
1876+ for (ib = 0; ib < nb; ib++)
1877+ for (jb = 0; jb < nb; jb++)
1878+ for (ie = 0; ie < ne; ie++)
1879+ wght[ik][ib][jb][ie] = 0.0;
1880+
1881+ thr = 1.0e-10;
1882+
1883+ for (it = 0; it < 6*nk; it++) {
1884+
1885+ for (i4 = 0; i4 < 4; i4++)
1886+ for (ib = 0; ib < nb; ib++) {
1887+ ei1[i4][ib] = 0.0;
1888+ ej1[i4][ib] = 0.0;
1889+ }
1890+ for (i20 = 0; i20 < 20; i20++) {
1891+ for (i4 = 0; i4 < 4; i4++) {
1892+ for (ib = 0; ib < nb; ib++) {
1893+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
1894+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
1895+ }
1896+ }
1897+ }
1898+
1899+ for (ib = 0; ib < nb; ib++) {
1900+
1901+ for (i4 = 0; i4 < 4; i4++)
1902+ for (jb = 0; jb < nb; jb++)
1903+ for (ie = 0; ie < ne; ie++)
1904+ w1[ib][i4][jb][ie] = 0.0;
1905+
1906+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
1907+ eig_sort(4, e, indx);
1908+
1909+ if (e[0] <= 0.0 && 0.0 < e[1]) {
1910+
1911+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
1912+
1913+ if (v > thr) {
1914+ for (j4 = 0; j4 < 4; j4++) {
1915+ ei2[j4] = 0.0;
1916+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1917+ }
1918+ for (i4 = 0; i4 < 4; i4++)
1919+ for (j4 = 0; j4 < 4; j4++) {
1920+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1921+ for (jb = 0; jb < nb; jb++)
1922+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1923+ }
1924+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
1925+ for (i4 = 0; i4 < 4; i4++)
1926+ for (jb = 0; jb < nb; jb++)
1927+ for (j4 = 0; j4 < 4; j4++)
1928+ for (ie = 0; ie < ne; ie++)
1929+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
1930+ }
1931+ }
1932+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
1933+
1934+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
1935+
1936+ if (v > thr) {
1937+ for (j4 = 0; j4 < 4; j4++) {
1938+ ei2[j4] = 0.0;
1939+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1940+ }
1941+ for (i4 = 0; i4 < 4; i4++)
1942+ for (j4 = 0; j4 < 4; j4++) {
1943+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1944+ for (jb = 0; jb < nb; jb++)
1945+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1946+ }
1947+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
1948+ for (i4 = 0; i4 < 4; i4++)
1949+ for (jb = 0; jb < nb; jb++)
1950+ for (j4 = 0; j4 < 4; j4++)
1951+ for (ie = 0; ie < ne; ie++)
1952+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
1953+ }
1954+
1955+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
1956+
1957+ if (v > thr) {
1958+ for (j4 = 0; j4 < 4; j4++) {
1959+ ei2[j4] = 0.0;
1960+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1961+ }
1962+ for (i4 = 0; i4 < 4; i4++)
1963+ for (j4 = 0; j4 < 4; j4++) {
1964+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1965+ for (jb = 0; jb < nb; jb++)
1966+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1967+ }
1968+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
1969+ for (i4 = 0; i4 < 4; i4++)
1970+ for (jb = 0; jb < nb; jb++)
1971+ for (j4 = 0; j4 < 4; j4++)
1972+ for (ie = 0; ie < ne; ie++)
1973+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
1974+ }
1975+
1976+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
1977+
1978+ if (v > thr) {
1979+ for (j4 = 0; j4 < 4; j4++) {
1980+ ei2[j4] = 0.0;
1981+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
1982+ }
1983+ for (i4 = 0; i4 < 4; i4++)
1984+ for (j4 = 0; j4 < 4; j4++) {
1985+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
1986+ for (jb = 0; jb < nb; jb++)
1987+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
1988+ }
1989+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
1990+ for (i4 = 0; i4 < 4; i4++)
1991+ for (jb = 0; jb < nb; jb++)
1992+ for (j4 = 0; j4 < 4; j4++)
1993+ for (ie = 0; ie < ne; ie++)
1994+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
1995+ }
1996+ }
1997+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
1998+
1999+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
2000+
2001+ if (v > thr) {
2002+ for (j4 = 0; j4 < 4; j4++) {
2003+ ei2[j4] = 0.0;
2004+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
2005+ }
2006+ for (i4 = 0; i4 < 4; i4++)
2007+ for (j4 = 0; j4 < 4; j4++) {
2008+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
2009+ for (jb = 0; jb < nb; jb++)
2010+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
2011+ }
2012+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
2013+ for (i4 = 0; i4 < 4; i4++)
2014+ for (jb = 0; jb < nb; jb++)
2015+ for (j4 = 0; j4 < 4; j4++)
2016+ for (ie = 0; ie < ne; ie++)
2017+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
2018+ }
2019+
2020+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
2021+
2022+ if (v > thr) {
2023+ for (j4 = 0; j4 < 4; j4++) {
2024+ ei2[j4] = 0.0;
2025+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
2026+ }
2027+ for (i4 = 0; i4 < 4; i4++)
2028+ for (j4 = 0; j4 < 4; j4++) {
2029+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
2030+ for (jb = 0; jb < nb; jb++)
2031+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
2032+ }
2033+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
2034+ for (i4 = 0; i4 < 4; i4++)
2035+ for (jb = 0; jb < nb; jb++)
2036+ for (j4 = 0; j4 < 4; j4++)
2037+ for (ie = 0; ie < ne; ie++)
2038+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
2039+ }
2040+
2041+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
2042+
2043+ if (v > thr) {
2044+ for (j4 = 0; j4 < 4; j4++) {
2045+ ei2[j4] = 0.0;
2046+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
2047+ }
2048+ for (i4 = 0; i4 < 4; i4++)
2049+ for (j4 = 0; j4 < 4; j4++) {
2050+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
2051+ for (jb = 0; jb < nb; jb++)
2052+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
2053+ }
2054+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
2055+ for (i4 = 0; i4 < 4; i4++)
2056+ for (jb = 0; jb < nb; jb++)
2057+ for (j4 = 0; j4 < 4; j4++)
2058+ for (ie = 0; ie < ne; ie++)
2059+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
2060+ }
2061+ }
2062+ else if (e[3] <= 0.0) {
2063+ for (i4 = 0; i4 < 4; i4++) {
2064+ ei2[i4] = ei1[i4][ib];
2065+ for (jb = 0; jb < nb; jb++)
2066+ ej2[jb][i4] = ej1[i4][jb];
2067+ }
2068+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
2069+ for (i4 = 0; i4 < 4; i4++)
2070+ for (jb = 0; jb < nb; jb++)
2071+ for (ie = 0; ie < ne; ie++)
2072+ w1[ib][i4][jb][ie] += w2[jb][i4][ie];
2073+ }
2074+ else {
2075+ continue;
2076+ }
2077+ }
2078+ for (i20 = 0; i20 < 20; i20++)
2079+ for (ib = 0; ib < nb; ib++)
2080+ for (i4 = 0; i4 < 4; i4++)
2081+ for (jb = 0; jb < nb; jb++)
2082+ for (ie = 0; ie < ne; ie++)
2083+ wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][i4] * w1[ib][i4][jb][ie];
2084+ }
2085+ for (ik = 0; ik < nk; ik++)
2086+ for (ib = 0; ib < nb; ib++)
2087+ for (jb = 0; jb < nb; jb++)
2088+ for (ie = 0; ie < ne; ie++)
2089+ wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
2090+
2091+ free(ikv[0]);
2092+ free(ikv);
2093+ free(ei1[0]);
2094+ free(ei1);
2095+ free(ej1[0]);
2096+ free(ej1);
2097+ free(ej2[0]);
2098+ free(ej2);
2099+ free(w1[0][0][0]);
2100+ free(w1[0][0]);
2101+ free(w1[0]);
2102+ free(w1);
2103+ free(w2[0][0]);
2104+ free(w2[0]);
2105+ free(w2);
2106+ /*
2107+ Convert weight to python list object
2108+ */
2109+ wght_po = PyList_New(nk * nb * nb * ne);
2110+ ierr = 0;
2111+ for (ik = 0; ik < nk; ik++)
2112+ for (ib = 0; ib < nb; ib++)
2113+ for (jb = 0; jb < nb; jb++)
2114+ for (ie = 0; ie < ne; ie++)
2115+ ierr = PyList_SetItem(wght_po, ik * nb * nb * ne + ib * nb * ne + jb * ne + ie,
2116+ PyFloat_FromDouble(wght[ik][ib][jb][ie]));
2117+ if (ierr != 0)printf("Error in PyList_SetItem\n");
2118+
2119+ free(eig1[0]);
2120+ free(eig1);
2121+ free(eig2[0]);
2122+ free(eig2);
2123+ free(wght[0][0][0]);
2124+ free(wght[0][0]);
2125+ free(wght[0]);
2126+ free(wght);
2127+ free(e0);
2128+
2129+ return wght_po;
2130+}
2131+/*
2132+Main SUBROUTINE for occupation : Theta(EF - E1)
2133+*/
2134+void occ_main(
2135+ int nk,
2136+ int nb,
2137+ double** eig,
2138+ int **ikv,
2139+ double wlsm[20][4],
2140+ double **ei1,
2141+ double **w1,
2142+ double **wght
2143+) {
2144+ int it, ik, ib, ii, jj, indx[4];
2145+ double e[4], v, tsmall[4][4];
2146+
2147+ for (ik = 0; ik < nk; ik++)
2148+ for (ib = 0; ib < nb; ib++)
2149+ wght[ik][ib] = 0.0;
2150+
2151+ for (it = 0; it < 6 * nk; it++) {
2152+
2153+ for (ii = 0; ii < 4; ii++)
2154+ for (ib = 0; ib < nb; ib++)
2155+ ei1[ii][ib] = 0.0;
2156+ for (jj = 0; jj < 20; jj++)
2157+ for (ii = 0; ii < 4; ii++)
2158+ for (ib = 0; ib < nb; ib++)
2159+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
2160+
2161+ for (ib = 0; ib < nb; ib++) {
2162+
2163+ for (ii = 0; ii < 4; ii++)
2164+ w1[ib][ii] = 0.0;
2165+
2166+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
2167+ eig_sort(4, e, indx);
2168+
2169+ if (e[0] <= 0.0 && 0.0 < e[1]) {
2170+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
2171+ for (ii = 0; ii < 4; ii++)
2172+ for (jj = 0; jj < 4; jj++)
2173+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2174+ }
2175+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
2176+
2177+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
2178+ for (ii = 0; ii < 4; ii++)
2179+ for (jj = 0; jj < 4; jj++)
2180+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2181+
2182+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
2183+ for (ii = 0; ii < 4; ii++)
2184+ for (jj = 0; jj < 4; jj++)
2185+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2186+
2187+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
2188+ for (ii = 0; ii < 4; ii++)
2189+ for (jj = 0; jj < 4; jj++)
2190+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2191+ }
2192+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
2193+
2194+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
2195+ for (ii = 0; ii < 4; ii++)
2196+ for (jj = 0; jj < 4; jj++)
2197+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2198+
2199+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
2200+ for (ii = 0; ii < 4; ii++)
2201+ for (jj = 0; jj < 4; jj++)
2202+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2203+
2204+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
2205+ for (ii = 0; ii < 4; ii++)
2206+ for (jj = 0; jj < 4; jj++)
2207+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
2208+ }
2209+ else if (e[3] <= 0.0) {
2210+ for (ii = 0; ii < 4; ii++)
2211+ w1[ib][ii] += 0.25;
2212+ }
2213+ else {
2214+ continue;
2215+ }
2216+ }
2217+ for (ii = 0; ii < 20; ii++)
2218+ for (ib = 0; ib < nb; ib++)
2219+ for (jj = 0; jj < 4; jj++) {
2220+ wght[ikv[it][ii]][ib] += wlsm[ii][jj] * w1[ib][jj];
2221+ }
2222+ }
2223+ for (ik = 0; ik < nk; ik++)
2224+ for (ib = 0; ib < nb; ib++)
2225+ wght[ik][ib] /= (6.0 * (double)nk);
2226+}
2227+/*
2228+
2229+*/
2230+static PyObject* occ(PyObject* self, PyObject* args)
2231+{
2232+ int ik, ib, ii, nk, nb, ng[3], ierr, ** ikv;
2233+ double** eig, ** wght, bvec[3][3], wlsm[20][4], **ei1, **w1;
2234+ PyObject* eig_po, * wght_po;
2235+ /*
2236+ Read input from python object
2237+ */
2238+ if (!PyArg_ParseTuple(args, "iiiiidddddddddO",
2239+ &ng[0], &ng[1], &ng[2], &nk, &nb,
2240+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
2241+ &eig_po))
2242+ return NULL;
2243+ /*
2244+ convert python list object to array
2245+ */
2246+ eig = (double**)malloc(nk * sizeof(double*));
2247+ eig[0] = (double*)malloc(nk * nb * sizeof(double));
2248+ wght = (double**)malloc(nk * sizeof(double*));
2249+ wght[0] = (double*)malloc(nk * nb * sizeof(double));
2250+ for (ik = 0; ik < nk; ik++) {
2251+ eig[ik] = eig[0] + ik * nb;
2252+ wght[ik] = wght[0] + ik * nb;
2253+ }
2254+
2255+ for (ik = 0; ik < nk; ik++)
2256+ for (ib = 0; ib < nb; ib++)
2257+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
2258+
2259+ /*
2260+ Start main calculation
2261+ */
2262+ ikv = (int**)malloc(6 * nk * sizeof(int*));
2263+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
2264+ for (ik = 0; ik < 6 * nk; ik++) {
2265+ ikv[ik] = ikv[0] + ik * 20;
2266+ }
2267+ ei1 = (double**)malloc(4 * sizeof(double*));
2268+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
2269+ for (ii = 0; ii < 4; ii++) {
2270+ ei1[ii] = ei1[0] + ii * nb;
2271+ }
2272+
2273+ w1 = (double**)malloc(nb * sizeof(double*));
2274+ w1[0] = (double*)malloc(nb * 4 * sizeof(double));
2275+ for (ib = 0; ib < nb; ib++) {
2276+ w1[ib] = w1[0] + ib * 4;
2277+ }
2278+
2279+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
2280+
2281+ occ_main(nk, nb, eig, ikv, wlsm, ei1, w1, wght);
2282+
2283+ free(ikv[0]);
2284+ free(ikv);
2285+ free(ei1[0]);
2286+ free(ei1);
2287+ free(w1[0]);
2288+ free(w1);
2289+ /*
2290+ Convert weight to python list object
2291+ */
2292+ wght_po = PyList_New(nk * nb);
2293+ ierr = 0;
2294+ for (ik = 0; ik < nk; ik++)
2295+ for (ib = 0; ib < nb; ib++)
2296+ ierr = PyList_SetItem(wght_po, ik * nb + ib, PyFloat_FromDouble(wght[ik][ib]));
2297+ if (ierr != 0)printf("Error in PyList_SetItem\n");
2298+
2299+ free(eig[0]);
2300+ free(eig);
2301+ free(wght[0]);
2302+ free(wght);
2303+ return wght_po;
2304+}
2305+/*
2306+Calculate Fermi energy
2307+*/
2308+static PyObject* fermieng(PyObject* self, PyObject* args)
2309+{
2310+ int maxiter, ik, ib, iteration, nk, nb, ng[3],
2311+ **ikv, ii, ierr;
2312+ double eps, elw, eup, **eig2, sumkmid, bvec[3][3],
2313+ **eig, **wght, wlsm[20][4], **ei1, **w1, nelec, ef;
2314+ PyObject* eig_po, * wght_po;
2315+ /*
2316+ Read input from python object
2317+ */
2318+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOd",
2319+ &ng[0], &ng[1], &ng[2], &nk, &nb,
2320+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
2321+ &eig_po, &nelec))
2322+ return NULL;
2323+ /*
2324+ convert python list object to array
2325+ */
2326+ eig = (double**)malloc(nk * sizeof(double*));
2327+ eig[0] = (double*)malloc(nk * nb * sizeof(double));
2328+ wght = (double**)malloc(nk * sizeof(double*));
2329+ wght[0] = (double*)malloc(nk * nb * sizeof(double));
2330+ for (ik = 0; ik < nk; ik++) {
2331+ eig[ik] = eig[0] + ik * nb;
2332+ wght[ik] = wght[0] + ik * nb;
2333+ }
2334+
2335+ for (ik = 0; ik < nk; ik++)
2336+ for (ib = 0; ib < nb; ib++)
2337+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
2338+
2339+ /*
2340+ Start main calculation
2341+ */
2342+ ikv = (int**)malloc(6 * nk * sizeof(int*));
2343+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
2344+ for (ik = 0; ik < 6 * nk; ik++) {
2345+ ikv[ik] = ikv[0] + ik * 20;
2346+ }
2347+ ei1 = (double**)malloc(4 * sizeof(double*));
2348+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
2349+ for (ii = 0; ii < 4; ii++) {
2350+ ei1[ii] = ei1[0] + ii * nb;
2351+ }
2352+
2353+ w1 = (double**)malloc(nb * sizeof(double*));
2354+ w1[0] = (double*)malloc(nb * 4 * sizeof(double));
2355+ for (ib = 0; ib < nb; ib++) {
2356+ w1[ib] = w1[0] + ib * 4;
2357+ }
2358+
2359+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
2360+
2361+ eig2 = (double**)malloc(nk * sizeof(double*));
2362+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
2363+ for (ik = 0; ik < nk; ik++) {
2364+ eig2[ik] = eig2[0] + ik * nb;
2365+ }
2366+
2367+ maxiter = 300;
2368+ eps = 1.0e-10;
2369+
2370+ elw = eig[0][0];
2371+ eup = eig[0][0];
2372+ for (ik = 0; ik < nk; ik++) {
2373+ for (ib = 0; ib < nb; ib++) {
2374+ if (elw > eig[ik][ib]) elw = eig[ik][ib];
2375+ if (eup < eig[ik][ib]) eup = eig[ik][ib];
2376+ }
2377+ }
2378+ /*
2379+ Bisection method
2380+ */
2381+ for (iteration = 0; iteration < maxiter; iteration++) {
2382+
2383+ ef = (eup + elw) * 0.5;
2384+ /*
2385+ Calc. # of electrons
2386+ */
2387+ for (ik = 0; ik < nk; ik++)
2388+ for (ib = 0; ib < nb; ib++)
2389+ eig2[ik][ib] = eig[ik][ib] - ef;
2390+ occ_main(nk, nb, eig2, ikv, wlsm, ei1, w1, wght);
2391+
2392+ sumkmid = 0.0;
2393+ for (ik = 0; ik < nk; ik++) {
2394+ for (ib = 0; ib < nb; ib++) {
2395+ sumkmid += wght[ik][ib];
2396+ }
2397+ }
2398+ /*
2399+ convergence check
2400+ */
2401+ if (fabs(sumkmid - nelec) < eps) {
2402+ break;
2403+ }
2404+ else if (sumkmid < nelec)
2405+ elw = ef;
2406+ else
2407+ eup = ef;
2408+ }
2409+ free(eig2);
2410+ free(ikv[0]);
2411+ free(ikv);
2412+ free(ei1[0]);
2413+ free(ei1);
2414+ free(w1[0]);
2415+ free(w1);
2416+ /*
2417+ Convert weight to python list object
2418+ */
2419+ wght_po = PyList_New(nk * nb+2);
2420+ ierr = 0;
2421+ for (ik = 0; ik < nk; ik++)
2422+ for (ib = 0; ib < nb; ib++)
2423+ ierr = PyList_SetItem(wght_po, ik * nb + ib, PyFloat_FromDouble(wght[ik][ib]));
2424+ ierr = PyList_SetItem(wght_po, nk*nb, PyFloat_FromDouble(ef));
2425+ ierr = PyList_SetItem(wght_po, nk * nb+1, PyLong_FromLong((long)iteration));
2426+ if (ierr != 0)printf("Error in PyList_SetItem\n");
2427+
2428+ free(eig[0]);
2429+ free(eig);
2430+ free(wght[0]);
2431+ free(wght);
2432+ return wght_po;
2433+}
2434+/*
2435+Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
2436+for 0<x<1, 0<y<1-x, 0<z<1-x-y
2437+*/
2438+/*
2439+1, Different each other
2440+*/
2441+static void libtetrabz_polcmplx_1234(
2442+ double g1,
2443+ double g2,
2444+ double g3,
2445+ double g4,
2446+ double* w
2447+) {
2448+ double w2, w3, w4;
2449+ /*
2450+ Real
2451+ */
2452+ w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
2453+ (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2454+ w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
2455+ w2 = w2 / (g2 - g1);
2456+ w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
2457+ (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2458+ w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
2459+ w3 = w3 / (g3 - g1);
2460+ w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
2461+ (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
2462+ w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
2463+ w4 = w4 / (g4 - g1);
2464+ w2 = (w2 - w3) / (g2 - g3);
2465+ w4 = (w4 - w3) / (g4 - g3);
2466+ w[0] = (w4 - w2) / (2.0 * (g4 - g2));
2467+ /*
2468+ Imaginal
2469+ */
2470+ w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
2471+ (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2472+ w2 = 4.0 * g2 - w2 / (g2 - g1);
2473+ w2 = w2 / (g2 - g1);
2474+ w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
2475+ (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2476+ w3 = 4.0 * g3 - w3 / (g3 - g1);
2477+ w3 = w3 / (g3 - g1);
2478+ w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
2479+ (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
2480+ w4 = 4.0 * g4 - w4 / (g4 - g1);
2481+ w4 = w4 / (g4 - g1);
2482+ w2 = (w2 - w3) / (g2 - g3);
2483+ w4 = (w4 - w3) / (g4 - g3);
2484+ w[1] = (w4 - w2) / (2.0 * (g4 - g2));
2485+}
2486+/*
2487+2, g4 = g1
2488+*/
2489+static void libtetrabz_polcmplx_1231(
2490+ double g1,
2491+ double g2,
2492+ double g3,
2493+ double* w
2494+) {
2495+ double w2, w3;
2496+ /*
2497+ Real
2498+ */
2499+ w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
2500+ + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2501+ w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
2502+ w2 = -g1 + w2 / (g2 - g1);
2503+ w2 = w2 / (g2 - g1);
2504+ w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
2505+ + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2506+ w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
2507+ w3 = -g1 + w3 / (g3 - g1);
2508+ w3 = w3 / (g3 - g1);
2509+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
2510+ /*
2511+ Imaginal
2512+ */
2513+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
2514+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2515+ w2 = 4.0 * g2 - w2 / (g2 - g1);
2516+ w2 = 1 + w2 / (g2 - g1);
2517+ w2 = w2 / (g2 - g1);
2518+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
2519+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2520+ w3 = 4.0 * g3 - w3 / (g3 - g1);
2521+ w3 = 1 + w3 / (g3 - g1);
2522+ w3 = w3 / (g3 - g1);
2523+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
2524+}
2525+/*
2526+3, g4 = g3
2527+*/
2528+static void libtetrabz_polcmplx_1233(
2529+ double g1,
2530+ double g2,
2531+ double g3,
2532+ double* w
2533+) {
2534+ double w2, w3;
2535+ /*
2536+ Real
2537+ */
2538+ w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
2539+ g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2540+ w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
2541+ w2 = w2 / (g2 - g1);
2542+ w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
2543+ g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2544+ w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
2545+ w3 = w3 / (g3 - g1);
2546+ w2 = (w3 - w2) / (g3 - g2);
2547+ w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
2548+ + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2549+ w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
2550+ w3 = 4.0 * g1 + w3 / (g3 - g1);
2551+ w3 = w3 / (g3 - g1);
2552+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
2553+ /*
2554+ Imaginal
2555+ */
2556+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
2557+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2558+ w2 = 4.0 * g2 - w2 / (g2 - g1);
2559+ w2 = w2 / (g2 - g1);
2560+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
2561+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2562+ w3 = 4.0 * g3 - w3 / (g3 - g1);
2563+ w3 = w3 / (g3 - g1);
2564+ w2 = (w3 - w2) / (g3 - g2);
2565+ w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3 * g3) * (atan(g3) - atan(g1))
2566+ + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
2567+ w3 = w3 / (g3 - g1) - 4.0 * g1;
2568+ w3 = w3 / (g3 - g1) - 2.0;
2569+ w3 = (2.0 * w3) / (g3 - g1);
2570+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
2571+}
2572+/*
2573+4, g4 = g1 and g3 = g2
2574+*/
2575+static void libtetrabz_polcmplx_1221(
2576+ double g1,
2577+ double g2,
2578+ double* w
2579+) {
2580+ /*
2581+ Real
2582+ */
2583+ w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
2584+ + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2585+ w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
2586+ w[0] = 3.0 * g1 + w[0] / (g2 - g1);
2587+ w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
2588+ w[0] = w[0] / (2.0 * (g2 - g1));
2589+ /*
2590+ Imaginal
2591+ */
2592+ w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
2593+ + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
2594+ w[1] = -4.0 * g1 + w[1] / (g2 - g1);
2595+ w[1] = -3.0 + w[1] / (g2 - g1);
2596+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
2597+}
2598+/*
2599+5, g4 = g3 = g2
2600+*/
2601+static void libtetrabz_polcmplx_1222(
2602+ double g1,
2603+ double g2,
2604+ double* w
2605+) {
2606+ /*
2607+ Real
2608+ */
2609+ w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
2610+ + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2611+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
2612+ w[0] = g1 - w[0] / (g2 - g1);
2613+ w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
2614+ w[0] = w[0] / (2.0 * (g2 - g1));
2615+ /*
2616+ Imaginal
2617+ */
2618+ w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
2619+ + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2620+ w[1] = 4.0 * g1 + w[1] / (g2 - g1);
2621+ w[1] = 1.0 + w[1] / (g2 - g1);
2622+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
2623+}
2624+/*
2625+6, g4 = g3 = g1
2626+*/
2627+static void libtetrabz_polcmplx_1211(
2628+ double g1,
2629+ double g2,
2630+ double* w
2631+) {
2632+ /*
2633+ Real
2634+ */
2635+ w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
2636+ + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2637+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
2638+ w[0] = -5.0 * g1 + w[0] / (g2 - g1);
2639+ w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
2640+ w[0] = w[0] / (6.0 * (g2 - g1));
2641+ /*
2642+ Imaginal
2643+ */
2644+ w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
2645+ + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
2646+ w[1] = 4.0 * g2 + w[1] / (g2 - g1);
2647+ w[1] = 1.0 + w[1] / (g2 - g1);
2648+ w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
2649+}
2650+/*
2651+Tetrahedron method for delta(om - ep + e)
2652+*/
2653+static void libtetrabz_polcmplx3(
2654+ int ne,
2655+ double** e0,
2656+ double de[4],
2657+ double*** w1
2658+) {
2659+ int i4, ir, indx[4], ie;
2660+ double e[4], x[4], thr, w2[4][2], denom;
2661+
2662+ for (i4 = 0; i4 < 3; i4++) e[i4] = de[i4];
2663+ eig_sort(4, e, indx);
2664+
2665+ for (ie = 0; ie < ne; ie++) {
2666+ /*
2667+ I am not sure which one is better.
2668+ The former is more stable.
2669+ The latter is more accurate ?
2670+ */
2671+ for (i4 = 0; i4 < 4; i4++) {
2672+ denom = (de[i4] + e0[ie][0]) * (de[i4] + e0[ie][0]) + e0[ie][1] * e0[ie][1];
2673+ w1[i4][ie][0] = 0.25 * (de[i4] + e0[ie][0]) / denom;
2674+ w1[i4][ie][1] = 0.25 * (-e0[ie][1]) / denom;
2675+ }
2676+ continue;
2677+
2678+ for (i4 = 0; i4 < 4; i4++)
2679+ x[i4] = (e[i4] + e0[ie][0]) / e0[ie][1];
2680+ /* thr = maxval(de(1:4)) * 1d-3*/
2681+ thr = fabs(x[0]);
2682+ for (i4 = 0; i4 < 4; i4++)
2683+ if (thr < fabs(x[i4])) thr = fabs(x[i4]);
2684+ thr = fmax(1.0e-3, thr * 1.0e-2);
2685+
2686+ if (fabs(x[3] - x[2]) < thr) {
2687+ if (fabs(x[3] - x[1]) < thr) {
2688+ if (fabs(x[3] - x[0]) < thr) {
2689+ /*
2690+ e[3] = e[2] = e[1] = e[0]
2691+ */
2692+ w2[3][0] = 0.25 * x[3] / (1.0 + x[3] * x[3]);
2693+ w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
2694+ for (ir = 0; ir < 2; ir++) {
2695+ w2[2][ir] = w2[3][ir];
2696+ w2[1][ir] = w2[3][ir];
2697+ w2[0][ir] = w2[3][ir];
2698+ }
2699+ }
2700+ else {
2701+ /*
2702+ e[3] = e[2] = e[1]
2703+ */
2704+ libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
2705+ for (ir = 0; ir < 2; ir++) {
2706+ w2[2][ir] = w2[3][ir];
2707+ w2[1][ir] = w2[3][ir];
2708+ }
2709+ libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
2710+ /*
2711+ # IF(ANY(w2(1:2,1:4) < 0.0)):
2712+ # WRITE(*,*) ie
2713+ # WRITE(*,'(100e15.5)') x[0:4]
2714+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
2715+ # STOP "weighting 4=3=2"
2716+ */
2717+ }
2718+ }
2719+ else if (fabs(x[1] - x[0]) < thr) {
2720+ /*
2721+ e[3] = e[2], e[1] = e[0]
2722+ */
2723+ libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
2724+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
2725+ libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
2726+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
2727+ /*
2728+ IF(ANY(w2(1:2,1:4) < 0.0)):
2729+ WRITE(*,*) ie
2730+ WRITE(*,'(100e15.5)') x[0:4]
2731+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2732+ STOP "weighting 4=3 2=1"
2733+ */
2734+ }
2735+ else {
2736+ /*
2737+ e[3] = e[2]
2738+ */
2739+ libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
2740+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
2741+ libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
2742+ libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
2743+ /*
2744+ IF(ANY(w2(1:2,1:4) < 0.0)):
2745+ WRITE(*,*) ie
2746+ WRITE(*,'(100e15.5)') x[0:4]
2747+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2748+ STOP "weighting 4=3"
2749+ */
2750+ }
2751+ }
2752+ else if (fabs(x[2] - x[1]) < thr) {
2753+ if (fabs(x[2] - x[0]) < thr) {
2754+ /*
2755+ e[2] = e[1] = e[0]
2756+ */
2757+ libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
2758+ libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
2759+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
2760+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[2][ir];
2761+ /*
2762+ IF(ANY(w2(1:2,1:4) < 0.0)):
2763+ WRITE(*,*) ie
2764+ WRITE(*,'(100e15.5)') x[0:4]
2765+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2766+ STOP "weighting 3=2=1"
2767+ */
2768+ }
2769+ else {
2770+ /*
2771+ e[2] = e[1]
2772+ */
2773+ libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
2774+ libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
2775+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
2776+ libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
2777+ /*
2778+ IF(ANY(w2(1:2,1:4) < 0.0)):
2779+ WRITE(*,*) ie
2780+ WRITE(*,'(100e15.5)') x[0:4]
2781+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2782+ STOP "weighting 3=2"
2783+ */
2784+ }
2785+ }
2786+ else if (fabs(x[1] - x[0]) < thr) {
2787+ /*
2788+ e[1] = e[0]
2789+ */
2790+ libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
2791+ libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
2792+ libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
2793+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
2794+ /*
2795+ IF(ANY(w2(1:2,1:4) < 0.0)):
2796+ WRITE(*,*) ie
2797+ WRITE(*,'(100e15.5)') x[0:4]
2798+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2799+ STOP "weighting 2=1"
2800+ */
2801+ }
2802+ else {
2803+ /*
2804+ Different each other.
2805+ */
2806+ libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2], w2[3]);
2807+ libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3], w2[2]);
2808+ libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3], w2[1]);
2809+ libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3], w2[0]);
2810+ /*
2811+ IF(ANY(w2(1:2,1:4) < 0.0)):
2812+ WRITE(*,*) ie
2813+ WRITE(*,'(100e15.5)') x[0:4]
2814+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
2815+ STOP "weighting"
2816+ */
2817+ }
2818+ for (i4 = 0; i4 < 4; i4++) {
2819+ w1[indx[i4]][ie][0] = w2[i4][0] / e0[ie][1];
2820+ w1[indx[i4]][ie][1] = w2[i4][1] / (-e0[ie][1]);
2821+ }
2822+ }
2823+}
2824+/*
2825+Tetrahedron method for theta( - E2)
2826+*/
2827+static void libtetrabz_polcmplx2(
2828+ int nb,
2829+ int ne,
2830+ double* *e0,
2831+ double* ei1,
2832+ double** ej1,
2833+ double**** w1
2834+) {
2835+ int ib, i4, j4, ir, ie, indx[4];
2836+ double e[4], tsmall[4][4], v, de[4], thr, *** w2;
2837+
2838+ w2 = (double***)malloc(4 * sizeof(double**));
2839+ w2[0] = (double**)malloc(4 * ne * sizeof(double*));
2840+ w2[0][0] = (double*)malloc(4 * ne * 2 * sizeof(double));
2841+ for (i4 = 0; i4 < 4; i4++) {
2842+ w2[i4] = w2[0] + i4 * ne;
2843+ for (ie = 0; ie < ne; ie++) {
2844+ w2[i4][ie] = w2[0][0] + i4 * ne * 2 + ie * 2;
2845+ }
2846+ }
2847+
2848+ thr = 1.0e-8;
2849+
2850+ for (ib = 0; ib < nb; ib++) {
2851+
2852+ for (i4 = 0; i4 < 4; i4++)
2853+ for (ie = 0; ie < ne; ie++)
2854+ for (ir = 0; ir < 2; ir++)
2855+ w1[ib][i4][ie][ir] = 0.0;
2856+
2857+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
2858+ eig_sort(4, e, indx);
2859+
2860+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
2861+
2862+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
2863+
2864+ if (v > thr) {
2865+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2866+ for (i4 = 0; i4 < 4; i4++)
2867+ for (j4 = 0; j4 < 4; j4++)
2868+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2869+ libtetrabz_polcmplx3(ne, e0, de, w2);
2870+ for (i4 = 0; i4 < 4; i4++)
2871+ for (j4 = 0; j4 < 4; j4++)
2872+ for (ie = 0; ie < ne; ie++)
2873+ for (ir = 0; ir < 2; ir++)
2874+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2875+ }
2876+ }
2877+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
2878+
2879+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
2880+
2881+ if (v > thr) {
2882+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2883+ for (i4 = 0; i4 < 4; i4++)
2884+ for (j4 = 0; j4 < 4; j4++)
2885+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2886+ libtetrabz_polcmplx3(ne, e0, de, w2);
2887+ for (i4 = 0; i4 < 4; i4++)
2888+ for (j4 = 0; j4 < 4; j4++)
2889+ for (ie = 0; ie < ne; ie++)
2890+ for (ir = 0; ir < 2; ir++)
2891+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2892+ }
2893+
2894+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
2895+
2896+ if (v > thr) {
2897+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2898+ for (i4 = 0; i4 < 4; i4++)
2899+ for (j4 = 0; j4 < 4; j4++)
2900+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2901+ libtetrabz_polcmplx3(ne, e0, de, w2);
2902+ for (i4 = 0; i4 < 4; i4++)
2903+ for (j4 = 0; j4 < 4; j4++)
2904+ for (ie = 0; ie < ne; ie++)
2905+ for (ir = 0; ir < 2; ir++)
2906+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2907+ }
2908+
2909+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
2910+
2911+ if (v > thr) {
2912+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2913+ for (i4 = 0; i4 < 4; i4++)
2914+ for (j4 = 0; j4 < 4; j4++)
2915+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2916+ libtetrabz_polcmplx3(ne, e0, de, w2);
2917+ for (i4 = 0; i4 < 4; i4++)
2918+ for (j4 = 0; j4 < 4; j4++)
2919+ for (ie = 0; ie < ne; ie++)
2920+ for (ir = 0; ir < 2; ir++)
2921+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2922+ }
2923+ }
2924+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
2925+
2926+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
2927+
2928+ if (v > thr) {
2929+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2930+ for (i4 = 0; i4 < 4; i4++)
2931+ for (j4 = 0; j4 < 4; j4++)
2932+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2933+ libtetrabz_polcmplx3(ne, e0, de, w2);
2934+ for (i4 = 0; i4 < 4; i4++)
2935+ for (j4 = 0; j4 < 4; j4++)
2936+ for (ie = 0; ie < ne; ie++)
2937+ for (ir = 0; ir < 2; ir++)
2938+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2939+ }
2940+
2941+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
2942+
2943+ if (v > thr) {
2944+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2945+ for (i4 = 0; i4 < 4; i4++)
2946+ for (j4 = 0; j4 < 4; j4++)
2947+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2948+ libtetrabz_polcmplx3(ne, e0, de, w2);
2949+ for (i4 = 0; i4 < 4; i4++)
2950+ for (j4 = 0; j4 < 4; j4++)
2951+ for (ie = 0; ie < ne; ie++)
2952+ for (ir = 0; ir < 2; ir++)
2953+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2954+ }
2955+
2956+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
2957+
2958+ if (v > thr) {
2959+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
2960+ for (i4 = 0; i4 < 4; i4++)
2961+ for (j4 = 0; j4 < 4; j4++)
2962+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
2963+ libtetrabz_polcmplx3(ne, e0, de, w2);
2964+ for (i4 = 0; i4 < 4; i4++)
2965+ for (j4 = 0; j4 < 4; j4++)
2966+ for (ie = 0; ie < ne; ie++)
2967+ for (ir = 0; ir < 2; ir++)
2968+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
2969+ }
2970+ }
2971+ else if (e[3] <= 0.0) {
2972+ for (i4 = 0; i4 < 4; i4++)
2973+ de[i4] = ej1[ib][i4] - ei1[i4];
2974+ libtetrabz_polcmplx3(ne, e0, de, w2);
2975+ for (i4 = 0; i4 < 4; i4++)
2976+ for (ie = 0; ie < ne; ie++)
2977+ for (ir = 0; ir < 2; ir++)
2978+ w1[ib][i4][ie][ir] += w2[i4][ie][ir];
2979+ }
2980+ }
2981+ free(w2[0][0]);
2982+ free(w2[0]);
2983+ free(w2);
2984+}
2985+/*
2986+Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
2987+*/
2988+static PyObject* polcmplx(PyObject* self, PyObject* args) {
2989+ int it, ik, ie, ib, i4, j4, ir, jb, **ikv, indx[4], i20, ierr, ng[3], nk, nb, ne;
2990+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], ***** w1, **** w2, v, tsmall[4][4], thr,
2991+ bvec[3][3], ** eig1, ** eig2, ** e0, ***** wght;
2992+ PyObject* eig1_po, * eig2_po, * e0_po, * wght_po;
2993+ /*
2994+ Read input from python object
2995+ */
2996+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOOO",
2997+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
2998+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
2999+ &eig1_po, &eig2_po, &e0_po))
3000+ return NULL;
3001+ /*
3002+ convert python list object to array
3003+ */
3004+ eig1 = (double**)malloc(nk * sizeof(double*));
3005+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
3006+ eig2 = (double**)malloc(nk * sizeof(double*));
3007+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
3008+ wght = (double*****)malloc(nk * sizeof(double****));
3009+ wght[0] = (double****)malloc(nk * nb * sizeof(double***));
3010+ wght[0][0] = (double***)malloc(nk * nb * nb * sizeof(double**));
3011+ wght[0][0][0] = (double**)malloc(nk * nb * nb * ne * sizeof(double*));
3012+ wght[0][0][0][0] = (double*)malloc(nk * nb * nb * ne * 2 * sizeof(double));
3013+ for (ik = 0; ik < nk; ik++) {
3014+ eig1[ik] = eig1[0] + ik * nb;
3015+ eig2[ik] = eig2[0] + ik * nb;
3016+ wght[ik] = wght[0] + ik * nb;
3017+ for (ib = 0; ib < nb; ib++) {
3018+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
3019+ for (jb = 0; jb < nb; jb++) {
3020+ wght[ik][ib][jb] = wght[0][0][0]
3021+ + ik * nb * nb * ne + ib * nb * ne + jb * ne;
3022+ for (ie = 0; ie < ne; ie++) {
3023+ wght[ik][ib][jb][ie] = wght[0][0][0][0]
3024+ + ik * nb * nb * ne * 2 + ib * nb * ne * 2 + jb * ne * 2 + ie * 2;
3025+ }
3026+ }
3027+ }
3028+ }
3029+ e0 = (double**)malloc(ne * sizeof(double*));
3030+ e0[0] = (double*)malloc(ne * 2 * sizeof(double));
3031+ for (ie = 0; ie < ne; ie++)
3032+ e0[ie] = e0[0] + ie * 2;
3033+
3034+ for (ik = 0; ik < nk; ik++)
3035+ for (ib = 0; ib < nb; ib++) {
3036+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
3037+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
3038+ }
3039+ for (ie = 0; ie < ne; ie++)
3040+ for (ir = 0; ir < 2; ir++)
3041+ e0[ie][ir] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie * 2 + ir));
3042+ /*
3043+ Start main calculation
3044+ */
3045+ ikv = (int**)malloc(6 * nk * sizeof(int*));
3046+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
3047+ for (it = 0; it < 6 * nk; it++) {
3048+ ikv[it] = ikv[0] + it * 20;
3049+ }
3050+
3051+ ei1 = (double**)malloc(4 * sizeof(double*));
3052+ ej1 = (double**)malloc(4 * sizeof(double*));
3053+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
3054+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
3055+ for (i4 = 0; i4 < 4; i4++) {
3056+ ei1[i4] = ei1[0] + i4 * nb;
3057+ ej1[i4] = ej1[0] + i4 * nb;
3058+ }
3059+
3060+ ej2 = (double**)malloc(nb * sizeof(double*));
3061+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
3062+ for (ib = 0; ib < nb; ib++)
3063+ ej2[ib] = ej2[0] + ib * 4;
3064+
3065+ w1 = (double*****)malloc(nb * sizeof(double****));
3066+ w1[0] = (double****)malloc(nb * 4 * sizeof(double***));
3067+ w1[0][0] = (double***)malloc(nb * 4 * nb * sizeof(double**));
3068+ w1[0][0][0] = (double**)malloc(nb * 4 * nb * ne * sizeof(double*));
3069+ w1[0][0][0][0] = (double*)malloc(nb * 4 * nb * ne * 2 * sizeof(double));
3070+ for (ib = 0; ib < nb; ib++) {
3071+ w1[ib] = w1[0] + ib * 4;
3072+ for (i4 = 0; i4 < 4; i4++) {
3073+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
3074+ for (jb = 0; jb < nb; jb++) {
3075+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
3076+ for (ie = 0; ie < ne; ie++) {
3077+ w1[ib][i4][jb][ie] = w1[0][0][0][0] + ib * 4 * nb * ne * 2 + i4 * nb * ne * 2 + jb * ne * 2 + ie * 2;
3078+ }
3079+ }
3080+ }
3081+ }
3082+
3083+ w2 = (double****)malloc(nb * sizeof(double***));
3084+ w2[0] = (double***)malloc(nb * 4 * sizeof(double**));
3085+ w2[0][0] = (double**)malloc(nb * 4 * ne * sizeof(double*));
3086+ w2[0][0][0] = (double*)malloc(nb * 4 * ne * 2 * sizeof(double));
3087+ for (ib = 0; ib < nb; ib++) {
3088+ w2[ib] = w2[0] + ib * 4;
3089+ for (i4 = 0; i4 < 4; i4++) {
3090+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
3091+ for (ie = 0; ie < ne; ie++) {
3092+ w2[ib][i4][ie] = w2[0][0][0] + ib * 4 * ne * 2 + i4 * ne * 2 + ie * 2;
3093+ }
3094+ }
3095+ }
3096+
3097+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
3098+
3099+ for (ik = 0; ik < nk; ik++)
3100+ for (ib = 0; ib < nb; ib++)
3101+ for (jb = 0; jb < nb; jb++)
3102+ for (ie = 0; ie < ne; ie++)
3103+ for (ir = 0; ir < 2; ir++)
3104+ wght[ik][ib][jb][ie][ir] = 0.0;
3105+
3106+ thr = 1.0e-8;
3107+
3108+ for (it = 0; it < 6 * nk; it++) {
3109+
3110+ for (i4 = 0; i4 < 4; i4++)
3111+ for (ib = 0; ib < nb; ib++) {
3112+ ei1[i4][ib] = 0.0;
3113+ ej1[i4][ib] = 0.0;
3114+ }
3115+ for (i20 = 0; i20 < 20; i20++) {
3116+ for (i4 = 0; i4 < 4; i4++) {
3117+ for (ib = 0; ib < nb; ib++) {
3118+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
3119+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
3120+ }
3121+ }
3122+ }
3123+
3124+ for (ib = 0; ib < nb; ib++) {
3125+
3126+ for (i4 = 0; i4 < 4; i4++)
3127+ for (jb = 0; jb < nb; jb++)
3128+ for (ie = 0; ie < ne; ie++)
3129+ for (ir = 0; ir < 2; ir++)
3130+ w1[ib][i4][jb][ie][ir] = 0.0;
3131+
3132+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
3133+ eig_sort(4, e, indx);
3134+
3135+ if (e[0] <= 0.0 && 0.0 < e[1]) {
3136+
3137+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
3138+
3139+ if (v > thr) {
3140+ for (j4 = 0; j4 < 4; j4++) {
3141+ ei2[j4] = 0.0;
3142+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3143+ }
3144+ for (i4 = 0; i4 < 4; i4++)
3145+ for (j4 = 0; j4 < 4; j4++) {
3146+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3147+ for (jb = 0; jb < nb; jb++)
3148+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3149+ }
3150+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3151+ for (i4 = 0; i4 < 4; i4++)
3152+ for (jb = 0; jb < nb; jb++)
3153+ for (j4 = 0; j4 < 4; j4++)
3154+ for (ie = 0; ie < ne; ie++)
3155+ for (ir = 0; ir < 2; ir++)
3156+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3157+ }
3158+ }
3159+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
3160+
3161+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
3162+
3163+ if (v > thr) {
3164+ for (j4 = 0; j4 < 4; j4++) {
3165+ ei2[j4] = 0.0;
3166+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3167+ }
3168+ for (i4 = 0; i4 < 4; i4++)
3169+ for (j4 = 0; j4 < 4; j4++) {
3170+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3171+ for (jb = 0; jb < nb; jb++)
3172+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3173+ }
3174+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3175+ for (i4 = 0; i4 < 4; i4++)
3176+ for (jb = 0; jb < nb; jb++)
3177+ for (j4 = 0; j4 < 4; j4++)
3178+ for (ie = 0; ie < ne; ie++)
3179+ for (ir = 0; ir < 2; ir++)
3180+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3181+ }
3182+
3183+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
3184+
3185+ if (v > thr) {
3186+ for (j4 = 0; j4 < 4; j4++) {
3187+ ei2[j4] = 0.0;
3188+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3189+ }
3190+ for (i4 = 0; i4 < 4; i4++)
3191+ for (j4 = 0; j4 < 4; j4++) {
3192+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3193+ for (jb = 0; jb < nb; jb++)
3194+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3195+ }
3196+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3197+ for (i4 = 0; i4 < 4; i4++)
3198+ for (jb = 0; jb < nb; jb++)
3199+ for (j4 = 0; j4 < 4; j4++)
3200+ for (ie = 0; ie < ne; ie++)
3201+ for (ir = 0; ir < 2; ir++)
3202+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3203+ }
3204+
3205+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
3206+
3207+ if (v > thr) {
3208+ for (j4 = 0; j4 < 4; j4++) {
3209+ ei2[j4] = 0.0;
3210+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3211+ }
3212+ for (i4 = 0; i4 < 4; i4++)
3213+ for (j4 = 0; j4 < 4; j4++) {
3214+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3215+ for (jb = 0; jb < nb; jb++)
3216+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3217+ }
3218+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3219+ for (i4 = 0; i4 < 4; i4++)
3220+ for (jb = 0; jb < nb; jb++)
3221+ for (j4 = 0; j4 < 4; j4++)
3222+ for (ie = 0; ie < ne; ie++)
3223+ for (ir = 0; ir < 2; ir++)
3224+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3225+ }
3226+ }
3227+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
3228+
3229+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
3230+
3231+ if (v > thr) {
3232+ for (j4 = 0; j4 < 4; j4++) {
3233+ ei2[j4] = 0.0;
3234+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3235+ }
3236+ for (i4 = 0; i4 < 4; i4++)
3237+ for (j4 = 0; j4 < 4; j4++) {
3238+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3239+ for (jb = 0; jb < nb; jb++)
3240+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3241+ }
3242+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3243+ for (i4 = 0; i4 < 4; i4++)
3244+ for (jb = 0; jb < nb; jb++)
3245+ for (j4 = 0; j4 < 4; j4++)
3246+ for (ie = 0; ie < ne; ie++)
3247+ for (ir = 0; ir < 2; ir++)
3248+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3249+ }
3250+
3251+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
3252+
3253+ if (v > thr) {
3254+ for (j4 = 0; j4 < 4; j4++) {
3255+ ei2[j4] = 0.0;
3256+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3257+ }
3258+ for (i4 = 0; i4 < 4; i4++)
3259+ for (j4 = 0; j4 < 4; j4++) {
3260+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3261+ for (jb = 0; jb < nb; jb++)
3262+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3263+ }
3264+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3265+ for (i4 = 0; i4 < 4; i4++)
3266+ for (jb = 0; jb < nb; jb++)
3267+ for (j4 = 0; j4 < 4; j4++)
3268+ for (ie = 0; ie < ne; ie++)
3269+ for (ir = 0; ir < 2; ir++)
3270+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3271+ }
3272+
3273+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
3274+
3275+ if (v > thr) {
3276+ for (j4 = 0; j4 < 4; j4++) {
3277+ ei2[j4] = 0.0;
3278+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3279+ }
3280+ for (i4 = 0; i4 < 4; i4++)
3281+ for (j4 = 0; j4 < 4; j4++) {
3282+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3283+ for (jb = 0; jb < nb; jb++)
3284+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3285+ }
3286+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3287+ for (i4 = 0; i4 < 4; i4++)
3288+ for (jb = 0; jb < nb; jb++)
3289+ for (j4 = 0; j4 < 4; j4++)
3290+ for (ie = 0; ie < ne; ie++)
3291+ for (ir = 0; ir < 2; ir++)
3292+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
3293+ }
3294+ }
3295+ else if (e[3] <= 0.0) {
3296+ for (i4 = 0; i4 < 4; i4++) {
3297+ ei2[i4] = ei1[i4][ib];
3298+ for (jb = 0; jb < nb; jb++)
3299+ ej2[jb][i4] = ej1[i4][jb];
3300+ }
3301+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
3302+ for (i4 = 0; i4 < 4; i4++)
3303+ for (jb = 0; jb < nb; jb++)
3304+ for (ie = 0; ie < ne; ie++)
3305+ for (ir = 0; ir < 2; ir++)
3306+ w1[ib][i4][jb][ie][ir] += w2[jb][i4][ie][ir];
3307+ }
3308+ else {
3309+ continue;
3310+ }
3311+ }
3312+ for (i20 = 0; i20 < 20; i20++)
3313+ for (ib = 0; ib < nb; ib++)
3314+ for (i4 = 0; i4 < 4; i4++)
3315+ for (jb = 0; jb < nb; jb++)
3316+ for (ie = 0; ie < ne; ie++)
3317+ for (ir = 0; ir < 2; ir++)
3318+ wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][i4] * w1[ib][i4][jb][ie][ir];
3319+ }
3320+ for (ik = 0; ik < nk; ik++)
3321+ for (ib = 0; ib < nb; ib++)
3322+ for (jb = 0; jb < nb; jb++)
3323+ for (ie = 0; ie < ne; ie++)
3324+ for (i4 = 0; i4 < 2; i4++)
3325+ wght[ik][ib][jb][ie][i4] /= (6.0 * (double) nk);
3326+
3327+ free(ikv[0]);
3328+ free(ikv);
3329+ free(ei1[0]);
3330+ free(ei1);
3331+ free(ej1[0]);
3332+ free(ej1);
3333+ free(ej2[0]);
3334+ free(ej2);
3335+ free(w1[0][0][0][0]);
3336+ free(w1[0][0][0]);
3337+ free(w1[0][0]);
3338+ free(w1[0]);
3339+ free(w1);
3340+ free(w2[0][0][0]);
3341+ free(w2[0][0]);
3342+ free(w2[0]);
3343+ free(w2);
3344+ /*
3345+ Convert weight to python list object
3346+ */
3347+ wght_po = PyList_New(nk * nb * nb * ne * 2);
3348+ ierr = 0;
3349+
3350+ for (ik = 0; ik < nk; ik++)
3351+ for (ib = 0; ib < nb; ib++)
3352+ for (jb = 0; jb < nb; jb++)
3353+ for (ie = 0; ie < ne; ie++)
3354+ for (ir = 0; ir < 2; ir++)
3355+ ierr = PyList_SetItem(wght_po,
3356+ ik * nb * nb * ne * 2 + ib * nb * ne * 2 + jb * ne * 2 + ie * 2 + ir,
3357+ PyFloat_FromDouble(wght[ik][ib][jb][ie][ir]));
3358+ if (ierr != 0)printf("Error in PyList_SetItem\n");
3359+
3360+ free(eig1[0]);
3361+ free(eig1);
3362+ free(eig2[0]);
3363+ free(eig2);
3364+ free(wght[0][0][0][0]);
3365+ free(wght[0][0][0]);
3366+ free(wght[0][0]);
3367+ free(wght[0]);
3368+ free(wght);
3369+ free(e0[0]);
3370+ free(e0);
3371+
3372+ return wght_po;
3373+}
3374+/*
3375+ Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
3376+ for 0<x<1, 0<y<1-x, 0<z<1-x-y
3377+*/
3378+/*
3379+1, Different each other
3380+*/
3381+static double libtetrabz_polstat_1234(
3382+ double g1,
3383+ double g2,
3384+ double g3,
3385+ double g4,
3386+ double lng1,
3387+ double lng2,
3388+ double lng3,
3389+ double lng4)
3390+{
3391+ double w2, w3, w4, w;
3392+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 / (g2 - g1);
3393+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 / (g3 - g1);
3394+ w4 = ((lng4 - lng1) / (g4 - g1) * g4 - 1.0) * g4 / (g4 - g1);
3395+ w2 = ((w2 - w3) * g2) / (g2 - g3);
3396+ w4 = ((w4 - w3) * g4) / (g4 - g3);
3397+ w = (w4 - w2) / (g4 - g2);
3398+ return w;
3399+}
3400+/*
3401+2, g4 = g1
3402+*/
3403+static double libtetrabz_polstat_1231(
3404+ double g1,
3405+ double g2,
3406+ double g3,
3407+ double lng1,
3408+ double lng2,
3409+ double lng3
3410+) {
3411+ double w2, w3, w;
3412+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 * g2 / (g2 - g1) - g1 / 2.0;
3413+ w2 = w2 / (g2 - g1);
3414+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 * g3 / (g3 - g1) - g1 / 2.0;
3415+ w3 = w3 / (g3 - g1);
3416+ w = (w3 - w2) / (g3 - g2);
3417+ return w;
3418+}
3419+/*
3420+# 3, g4 = g3
3421+*/
3422+static double libtetrabz_polstat_1233(
3423+ double g1,
3424+ double g2,
3425+ double g3,
3426+ double lng1,
3427+ double lng2,
3428+ double lng3
3429+) {
3430+ double w2, w3, w;
3431+ w2 = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
3432+ w2 = (g2 * w2) / (g2 - g1);
3433+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
3434+ w3 = (g3 * w3) / (g3 - g1);
3435+ w2 = (w3 - w2) / (g3 - g2);
3436+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
3437+ w3 = 1.0 - (2.0 * w3 * g1) / (g3 - g1);
3438+ w3 = w3 / (g3 - g1);
3439+ w = (g3 * w3 - g2 * w2) / (g3 - g2);
3440+ return w;
3441+}
3442+/*
3443+4, g4 = g1 and g3 = g2
3444+*/
3445+static double libtetrabz_polstat_1221(
3446+ double g1,
3447+ double g2,
3448+ double lng1,
3449+ double lng2
3450+) {
3451+ double w;
3452+ w = 1.0 - (lng2 - lng1) / (g2 - g1) * g1;
3453+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
3454+ w = -1.0 + (3.0 * g2 * w) / (g2 - g1);
3455+ w = w / (2.0 * (g2 - g1));
3456+ return w;
3457+}
3458+/*
3459+5, g4 = g3 = g2
3460+*/
3461+static double libtetrabz_polstat_1222(
3462+ double g1,
3463+ double g2,
3464+ double lng1,
3465+ double lng2
3466+) {
3467+ double w;
3468+ w = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
3469+ w = (2.0 * g1 * w) / (g2 - g1) - 1.0;
3470+ w = (3.0 * g1 * w) / (g2 - g1) + 1.0;
3471+ w = w / (2.0 * (g2 - g1));
3472+ return w;
3473+}
3474+/*
3475+6, g4 = g3 = g1
3476+*/
3477+static double libtetrabz_polstat_1211(
3478+ double g1,
3479+ double g2,
3480+ double lng1,
3481+ double lng2
3482+) {
3483+ double w;
3484+ w = -1.0 + (lng2 - lng1) / (g2 - g1) * g2;
3485+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
3486+ w = -1.0 + (3.0 * g2 * w) / (2.0 * (g2 - g1));
3487+ w = w / (3.0 * (g2 - g1));
3488+ return w;
3489+}
3490+/*
3491+Tetrahedron method for delta(om - ep + e)
3492+*/
3493+static void libtetrabz_polstat3(
3494+ double de[4],
3495+ double w1[4]
3496+)
3497+{
3498+ int i4, indx[4];
3499+ double thr, thr2, e[4], ln[4];
3500+
3501+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
3502+ eig_sort(4, e, indx);
3503+
3504+ thr = e[3] * 1.0e-3;
3505+ thr2 = 1.0e-8;
3506+
3507+ for(i4 =0; i4 <4; i4++){
3508+ if (e[i4] < thr2) {
3509+ if (i4 == 3) {
3510+ printf(" Nesting # \n");
3511+ }
3512+ ln[i4] = 0.0;
3513+ e[i4] = 0.0;
3514+ }
3515+ else{
3516+ ln[i4] = log(e[i4]);
3517+ }
3518+ }
3519+
3520+ if (fabs(e[3] - e[2]) < thr) {
3521+ if (fabs(e[3] - e[1]) < thr) {
3522+ if (fabs(e[3] - e[0]) < thr) {
3523+ /*
3524+ e[3] = e[2] = e[1] = e[0]
3525+ */
3526+ w1[indx[3]] = 0.25 / e[3];
3527+ w1[indx[2]] = w1[indx[3]];
3528+ w1[indx[1]] = w1[indx[3]];
3529+ w1[indx[0]] = w1[indx[3]];
3530+ }
3531+ else {
3532+ /*
3533+ e[3] = e[2] = e[1]
3534+ */
3535+ w1[indx[3]] = libtetrabz_polstat_1211(e[3], e[0], ln[3], ln[0]);
3536+ w1[indx[2]] = w1[indx[3]];
3537+ w1[indx[1]] = w1[indx[3]];
3538+ w1[indx[0]] = libtetrabz_polstat_1222(e[0], e[3], ln[0], ln[3]);
3539+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3540+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3541+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3542+ printf("weighting 4=3=2\n");
3543+ }
3544+ }
3545+ }
3546+ else if (fabs(e[1] - e[0]) < thr) {
3547+ /*
3548+ e[3] = e[2], e[1] = e[0]
3549+ */
3550+ w1[indx[3]] = libtetrabz_polstat_1221(e[3], e[1], ln[3], ln[1]);
3551+ w1[indx[2]] = w1[indx[3]];
3552+ w1[indx[1]] = libtetrabz_polstat_1221(e[1], e[3], ln[1], ln[3]);
3553+ w1[indx[0]] = w1[indx[1]];
3554+
3555+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3556+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3557+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3558+ printf("weighting 4=3 2=1\n");
3559+ }
3560+ }
3561+ else {
3562+ /*
3563+ e[3] = e[2]
3564+ */
3565+ w1[indx[3]] = libtetrabz_polstat_1231(e[3], e[0], e[1], ln[3], ln[0], ln[1]);
3566+ w1[indx[2]] = w1[indx[3]];
3567+ w1[indx[1]] = libtetrabz_polstat_1233(e[1], e[0], e[3], ln[1], ln[0], ln[3]);
3568+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[1], e[3], ln[0], ln[1], ln[3]);
3569+
3570+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3571+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3572+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3573+ printf("weighting 4=3\n");
3574+ }
3575+ }
3576+ }
3577+ else if (fabs(e[2] - e[1]) < thr) {
3578+ if (fabs(e[2] - e[0]) < thr) {
3579+ /*
3580+ e[2] = e[1] = e[0]
3581+ */
3582+ w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
3583+ w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
3584+ w1[indx[1]] = w1[indx[2]];
3585+ w1[indx[0]] = w1[indx[2]];
3586+
3587+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3588+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3589+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3590+ printf("weighting 3=2=1\n");
3591+ }
3592+ }
3593+ else {
3594+ /*
3595+ e[2] = e[1]
3596+ */
3597+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
3598+ w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
3599+ w1[indx[1]] = w1[indx[2]];
3600+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
3601+
3602+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3603+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3604+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3605+ printf("weighting 3=2\n");
3606+ }
3607+ }
3608+ }
3609+ else if (fabs(e[1] - e[0]) < thr) {
3610+ /*
3611+ e[1] = e[0]
3612+ */
3613+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
3614+ w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
3615+ w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
3616+ w1[indx[0]] = w1[indx[1]];
3617+
3618+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3619+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3620+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3621+ printf("weighting 2=1\n");
3622+ }
3623+ }
3624+ else {
3625+ /*
3626+ Different each other.
3627+ */
3628+ w1[indx[3]] = libtetrabz_polstat_1234(e[3], e[0], e[1], e[2], ln[3], ln[0], ln[1], ln[2]);
3629+ w1[indx[2]] = libtetrabz_polstat_1234(e[2], e[0], e[1], e[3], ln[2], ln[0], ln[1], ln[3]);
3630+ w1[indx[1]] = libtetrabz_polstat_1234(e[1], e[0], e[2], e[3], ln[1], ln[0], ln[2], ln[3]);
3631+ w1[indx[0]] = libtetrabz_polstat_1234(e[0], e[1], e[2], e[3], ln[0], ln[1], ln[2], ln[3]);
3632+
3633+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
3634+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
3635+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
3636+ printf("weighting\n");
3637+ }
3638+ }
3639+}
3640+/*
3641+Tetrahedron method for theta( - E2)
3642+*/
3643+static void libtetrabz_polstat2(
3644+ int nb,
3645+ double *ei1,
3646+ double **ej1,
3647+ double **w1
3648+) {
3649+ int i4, j4, ib, indx[4];
3650+ double v, thr, e[4], de[4], w2[4], tsmall[4][4];
3651+
3652+ thr = 1.0e-8;
3653+
3654+ for (ib = 0; ib < nb; ib++) {
3655+
3656+ for (i4 = 0; i4 < 4; i4++)
3657+ w1[ib][i4] = 0.0;
3658+
3659+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
3660+ eig_sort(4, e, indx);
3661+
3662+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
3663+
3664+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
3665+
3666+ if (v > thr) {
3667+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3668+ for (i4 = 0; i4 < 4; i4++)
3669+ for (j4 = 0; j4 < 4; j4++)
3670+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3671+ libtetrabz_polstat3(de, w2);
3672+ for (i4 = 0; i4 < 4; i4++)
3673+ for (j4 = 0; j4 < 4; j4++)
3674+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3675+ }
3676+ }
3677+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
3678+
3679+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
3680+
3681+ if (v > thr) {
3682+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3683+ for (i4 = 0; i4 < 4; i4++)
3684+ for (j4 = 0; j4 < 4; j4++)
3685+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3686+ libtetrabz_polstat3(de, w2);
3687+ for (i4 = 0; i4 < 4; i4++)
3688+ for (j4 = 0; j4 < 4; j4++)
3689+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3690+ }
3691+
3692+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
3693+
3694+ if (v > thr) {
3695+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3696+ for (i4 = 0; i4 < 4; i4++)
3697+ for (j4 = 0; j4 < 4; j4++)
3698+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3699+ libtetrabz_polstat3(de, w2);
3700+ for (i4 = 0; i4 < 4; i4++)
3701+ for (j4 = 0; j4 < 4; j4++)
3702+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3703+ }
3704+
3705+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
3706+
3707+ if (v > thr) {
3708+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3709+ for (i4 = 0; i4 < 4; i4++)
3710+ for (j4 = 0; j4 < 4; j4++)
3711+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3712+ libtetrabz_polstat3(de, w2);
3713+ for (i4 = 0; i4 < 4; i4++)
3714+ for (j4 = 0; j4 < 4; j4++)
3715+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3716+ }
3717+ }
3718+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
3719+
3720+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
3721+
3722+ if (v > thr) {
3723+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3724+ for (i4 = 0; i4 < 4; i4++)
3725+ for (j4 = 0; j4 < 4; j4++)
3726+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3727+ libtetrabz_polstat3(de, w2);
3728+ for (i4 = 0; i4 < 4; i4++)
3729+ for (j4 = 0; j4 < 4; j4++)
3730+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3731+ }
3732+
3733+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
3734+
3735+ if (v > thr) {
3736+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3737+ for (i4 = 0; i4 < 4; i4++)
3738+ for (j4 = 0; j4 < 4; j4++)
3739+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3740+ libtetrabz_polstat3(de, w2);
3741+ for (i4 = 0; i4 < 4; i4++)
3742+ for (j4 = 0; j4 < 4; j4++)
3743+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3744+ }
3745+
3746+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
3747+
3748+ if (v > thr) {
3749+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
3750+ for (i4 = 0; i4 < 4; i4++)
3751+ for (j4 = 0; j4 < 4; j4++)
3752+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
3753+ libtetrabz_polstat3(de, w2);
3754+ for (i4 = 0; i4 < 4; i4++)
3755+ for (j4 = 0; j4 < 4; j4++)
3756+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
3757+ }
3758+ }
3759+ else if (e[3] <= 0.0) {
3760+ for (i4 = 0; i4 < 4; i4++)
3761+ de[i4] = ej1[ib][i4] - ei1[i4];
3762+ libtetrabz_polstat3(de, w2);
3763+ for (i4 = 0; i4 < 4; i4++)
3764+ w1[ib][i4] += w2[i4];
3765+ }
3766+ }
3767+}
3768+/*
3769+Compute Static polarization function
3770+*/
3771+static PyObject* polstat(PyObject* self, PyObject* args)
3772+{
3773+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4], ierr, ng[3], nk, nb;
3774+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr,
3775+ bvec[3][3], ** eig1, ** eig2, *** wght;
3776+ PyObject* eig1_po, * eig2_po, * wght_po;
3777+ /*
3778+ Read input from python object
3779+ */
3780+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
3781+ &ng[0], &ng[1], &ng[2], &nk, &nb,
3782+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
3783+ &eig1_po, &eig2_po))
3784+ return NULL;
3785+ /*
3786+ convert python list object to array
3787+ */
3788+ eig1 = (double**)malloc(nk * sizeof(double*));
3789+ eig1[0] = (double*)malloc(nk * nb * sizeof(double));
3790+ eig2 = (double**)malloc(nk * sizeof(double*));
3791+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
3792+ wght = (double***)malloc(nk * sizeof(double**));
3793+ wght[0] = (double**)malloc(nk * nb * sizeof(double*));
3794+ wght[0][0] = (double*)malloc(nk * nb * nb * sizeof(double));
3795+ for (ik = 0; ik < nk; ik++) {
3796+ eig1[ik] = eig1[0] + ik * nb;
3797+ eig2[ik] = eig2[0] + ik * nb;
3798+ wght[ik] = wght[0] + ik * nb;
3799+ for (ib = 0; ib < nb; ib++) {
3800+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
3801+ }
3802+ }
3803+
3804+ for (ik = 0; ik < nk; ik++)
3805+ for (ib = 0; ib < nb; ib++) {
3806+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
3807+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
3808+ }
3809+ /*
3810+ Start main calculation
3811+ */
3812+
3813+
3814+ thr = 1.0e-10;
3815+
3816+ ikv = (int**)malloc(6 * nk * sizeof(int*));
3817+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
3818+ for (it = 0; it < 6 * nk; it++) {
3819+ ikv[it] = ikv[0] + it * 20;
3820+ }
3821+
3822+ ei1 = (double**)malloc(4 * sizeof(double*));
3823+ ej1 = (double**)malloc(4 * sizeof(double*));
3824+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
3825+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
3826+ for (i4 = 0; i4 < 4; i4++) {
3827+ ei1[i4] = ei1[0] + i4 * nb;
3828+ ej1[i4] = ej1[0] + i4 * nb;
3829+ }
3830+
3831+ ej2 = (double**)malloc(nb * sizeof(double*));
3832+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
3833+ for (ib = 0; ib < nb; ib++)
3834+ ej2[ib] = ej2[0] + ib * 4;
3835+
3836+ w1 = (double***)malloc(nb * sizeof(double**));
3837+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
3838+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
3839+ for (ib = 0; ib < nb; ib++) {
3840+ w1[ib] = w1[0] + ib * 4;
3841+ for (i4 = 0; i4 < 4; i4++) {
3842+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
3843+ }
3844+ }
3845+
3846+ w2 = (double**)malloc(nb * sizeof(double*));
3847+ w2[0] = (double*)malloc(nb * 4 * sizeof(double));
3848+ for (ib = 0; ib < nb; ib++) {
3849+ w2[ib] = w2[0] + ib * 4;
3850+ }
3851+
3852+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
3853+
3854+ for (ik = 0; ik < nk; ik++)
3855+ for (ib = 0; ib < nb; ib++)
3856+ for (jb = 0; jb < nb; jb++)
3857+ wght[ik][ib][jb] = 0.0;
3858+
3859+ for(it = 0; it < 6*nk; it++) {
3860+
3861+ for (i4 = 0; i4 < 4; i4++)
3862+ for (ib = 0; ib < nb; ib++) {
3863+ ei1[i4][ib] = 0.0;
3864+ ej1[i4][ib] = 0.0;
3865+ }
3866+ for (i20 = 0; i20 < 20; i20++) {
3867+ for (i4 = 0; i4 < 4; i4++) {
3868+ for (ib = 0; ib < nb; ib++) {
3869+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
3870+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
3871+ }
3872+ }
3873+ }
3874+
3875+ for (ib = 0; ib < nb; ib++) {
3876+
3877+ for (i4 = 0; i4 < 4; i4++)
3878+ for (jb = 0; jb < nb; jb++)
3879+ w1[ib][i4][jb] = 0.0;
3880+
3881+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
3882+ eig_sort(4, e, indx);
3883+
3884+ if (e[0] <= 0.0 && 0.0 < e[1]) {
3885+
3886+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
3887+
3888+ if (v > thr) {
3889+ for (j4 = 0; j4 < 4; j4++) {
3890+ ei2[j4] = 0.0;
3891+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3892+ }
3893+ for (i4 = 0; i4 < 4; i4++)
3894+ for (j4 = 0; j4 < 4; j4++) {
3895+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3896+ for (jb = 0; jb < nb; jb++)
3897+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3898+ }
3899+ libtetrabz_polstat2(nb, ei2, ej2, w2);
3900+ for (i4 = 0; i4 < 4; i4++)
3901+ for (jb = 0; jb < nb; jb++)
3902+ for (j4 = 0; j4 < 4; j4++)
3903+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
3904+ }
3905+ }
3906+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
3907+
3908+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
3909+
3910+ if (v > thr) {
3911+ for (j4 = 0; j4 < 4; j4++) {
3912+ ei2[j4] = 0.0;
3913+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3914+ }
3915+ for (i4 = 0; i4 < 4; i4++)
3916+ for (j4 = 0; j4 < 4; j4++) {
3917+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3918+ for (jb = 0; jb < nb; jb++)
3919+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3920+ }
3921+ libtetrabz_polstat2(nb, ei2, ej2, w2);
3922+ for (i4 = 0; i4 < 4; i4++)
3923+ for (jb = 0; jb < nb; jb++)
3924+ for (j4 = 0; j4 < 4; j4++)
3925+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
3926+ }
3927+
3928+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
3929+
3930+ if (v > thr) {
3931+ for (j4 = 0; j4 < 4; j4++) {
3932+ ei2[j4] = 0.0;
3933+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3934+ }
3935+ for (i4 = 0; i4 < 4; i4++)
3936+ for (j4 = 0; j4 < 4; j4++) {
3937+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3938+ for (jb = 0; jb < nb; jb++)
3939+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3940+ }
3941+ libtetrabz_polstat2(nb, ei2, ej2, w2);
3942+ for (i4 = 0; i4 < 4; i4++)
3943+ for (jb = 0; jb < nb; jb++)
3944+ for (j4 = 0; j4 < 4; j4++)
3945+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
3946+ }
3947+
3948+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
3949+
3950+ if (v > thr) {
3951+ for (j4 = 0; j4 < 4; j4++) {
3952+ ei2[j4] = 0.0;
3953+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3954+ }
3955+ for (i4 = 0; i4 < 4; i4++)
3956+ for (j4 = 0; j4 < 4; j4++) {
3957+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3958+ for (jb = 0; jb < nb; jb++)
3959+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3960+ }
3961+ libtetrabz_polstat2(nb, ei2, ej2, w2);
3962+ for (i4 = 0; i4 < 4; i4++)
3963+ for (jb = 0; jb < nb; jb++)
3964+ for (j4 = 0; j4 < 4; j4++)
3965+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
3966+ }
3967+ }
3968+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
3969+
3970+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
3971+
3972+ if (v > thr) {
3973+ for (j4 = 0; j4 < 4; j4++) {
3974+ ei2[j4] = 0.0;
3975+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3976+ }
3977+ for (i4 = 0; i4 < 4; i4++)
3978+ for (j4 = 0; j4 < 4; j4++) {
3979+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
3980+ for (jb = 0; jb < nb; jb++)
3981+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
3982+ }
3983+ libtetrabz_polstat2(nb, ei2, ej2, w2);
3984+ for (i4 = 0; i4 < 4; i4++)
3985+ for (jb = 0; jb < nb; jb++)
3986+ for (j4 = 0; j4 < 4; j4++)
3987+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
3988+ }
3989+
3990+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
3991+
3992+ if (v > thr) {
3993+ for (j4 = 0; j4 < 4; j4++) {
3994+ ei2[j4] = 0.0;
3995+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
3996+ }
3997+ for (i4 = 0; i4 < 4; i4++)
3998+ for (j4 = 0; j4 < 4; j4++) {
3999+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
4000+ for (jb = 0; jb < nb; jb++)
4001+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
4002+ }
4003+ libtetrabz_polstat2(nb, ei2, ej2, w2);
4004+ for (i4 = 0; i4 < 4; i4++)
4005+ for (jb = 0; jb < nb; jb++)
4006+ for (j4 = 0; j4 < 4; j4++)
4007+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
4008+ }
4009+
4010+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
4011+
4012+ if (v > thr) {
4013+ for (j4 = 0; j4 < 4; j4++) {
4014+ ei2[j4] = 0.0;
4015+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
4016+ }
4017+ for (i4 = 0; i4 < 4; i4++)
4018+ for (j4 = 0; j4 < 4; j4++) {
4019+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
4020+ for (jb = 0; jb < nb; jb++)
4021+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
4022+ }
4023+ libtetrabz_polstat2(nb, ei2, ej2, w2);
4024+ for (i4 = 0; i4 < 4; i4++)
4025+ for (jb = 0; jb < nb; jb++)
4026+ for (j4 = 0; j4 < 4; j4++)
4027+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
4028+ }
4029+ }
4030+ else if (e[3] <= 0.0) {
4031+ for (i4 = 0; i4 < 4; i4++) {
4032+ ei2[i4] = ei1[i4][ib];
4033+ for (jb = 0; jb < nb; jb++)
4034+ ej2[jb][i4] = ej1[i4][jb];
4035+ }
4036+ libtetrabz_polstat2(nb, ei2, ej2, w2);
4037+ for (i4 = 0; i4 < 4; i4++)
4038+ for (jb = 0; jb < nb; jb++)
4039+ w1[ib][i4][jb] += w2[jb][i4];
4040+ }
4041+ else {
4042+ continue;
4043+ }
4044+ }
4045+ for (i20 = 0; i20 < 20; i20++)
4046+ for (ib = 0; ib < nb; ib++)
4047+ for (i4 = 0; i4 < 4; i4++)
4048+ for (jb = 0; jb < nb; jb++)
4049+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
4050+ }
4051+ for (ik = 0; ik < nk; ik++)
4052+ for (ib = 0; ib < nb; ib++)
4053+ for (jb = 0; jb < nb; jb++)
4054+ wght[ik][ib][jb] /= (6.0 * (double) nk);
4055+
4056+ free(ikv[0]);
4057+ free(ikv);
4058+ free(ei1[0]);
4059+ free(ei1);
4060+ free(ej1[0]);
4061+ free(ej1);
4062+ free(ej2[0]);
4063+ free(ej2);
4064+ free(w1[0][0]);
4065+ free(w1[0]);
4066+ free(w1);
4067+ free(w2[0]);
4068+ free(w2);
4069+ /*
4070+ Convert weight to python list object
4071+ */
4072+ wght_po = PyList_New(nk * nb * nb);
4073+ ierr = 0;
4074+ for (ik = 0; ik < nk; ik++)
4075+ for (ib = 0; ib < nb; ib++)
4076+ for (jb = 0; jb < nb; jb++)
4077+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
4078+ if (ierr != 0)printf("Error in PyList_SetItem\n");
4079+
4080+ free(eig1[0]);
4081+ free(eig1);
4082+ free(eig2[0]);
4083+ free(eig2);
4084+ free(wght[0][0]);
4085+ free(wght[0]);
4086+ free(wght);
4087+
4088+ return wght_po;
4089+}
4090+static PyMethodDef LibTetraBZCMethods[] = {
4091+ {"dos", dos, METH_VARARGS, "Density of States"},
4092+ {"intdos", intdos, METH_VARARGS, "Integrated density of States"},
4093+ {"occ", occ, METH_VARARGS, "Occupations"},
4094+ {"fermieng", fermieng, METH_VARARGS, "Fermi energy"},
4095+ {"dbldelta", dbldelta, METH_VARARGS, "Double delta"},
4096+ {"dblstep", dblstep, METH_VARARGS, "Double step function"},
4097+ {"polstat", polstat, METH_VARARGS, "Static polarization function"},
4098+ {"fermigr", fermigr, METH_VARARGS, "Fermi's golden rule"},
4099+ {"polcmplx", polcmplx, METH_VARARGS, "Dynamical polarization function on complex frequency"},
4100+ {NULL, NULL, 0, NULL} /* Sentinel */
4101+};
4102+
4103+static struct PyModuleDef libtetrabzcmodule = {
4104+ PyModuleDef_HEAD_INIT,
4105+ "libtetrabzc", /* name of module */
4106+ NULL, /* module documentation, may be NULL */
4107+ -1, /* size of per-interpreter state of the module,
4108+ or -1 if the module keeps state in global variables. */
4109+ LibTetraBZCMethods
4110+};
4111+
4112+PyMODINIT_FUNC
4113+PyInit_libtetrabzc(void)
4114+{
4115+ return PyModule_Create(&libtetrabzcmodule);
4116+}
--- a/src/c/libtetrabz_common.c
+++ b/src/c/libtetrabz_common.c
@@ -21,6 +21,8 @@
2121 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 */
2323 #include <stdio.h>
24+#include <stdlib.h>
25+#include "libtetrabz.h"
2426
2527 void libtetrabz_initialize(
2628 int ng[3],
@@ -68,11 +70,11 @@ void libtetrabz_initialize(
6870 /*
6971 length of delta bvec
7072 */
71- for (i1 = 0; i1 < 4; i1++)
73+ for (i1 = 0; i1 < 4; i1++) {
7274 bnorm[i1] = 0.0;
73- for (i2 = 0; i2 < 3; i2++)
74- bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
75-
75+ for (i2 = 0; i2 < 3; i2++)
76+ bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
77+ }
7678 itype = 0;
7779 for (i1 = 1; i1 < 4; i1++)
7880 if (bnorm[i1] < bnorm[itype]) itype = i1;
--- a/src/c/libtetrabz_dbldelta.c
+++ b/src/c/libtetrabz_dbldelta.c
@@ -25,7 +25,7 @@
2525 #include <stdlib.h>
2626 #include <math.h>
2727
28-void libtetrabz_dbldelta2(
28+static void libtetrabz_dbldelta2(
2929 int nb,
3030 double **ej,
3131 double **w
--- a/src/c/libtetrabz_dblstep.c
+++ b/src/c/libtetrabz_dblstep.c
@@ -24,7 +24,7 @@
2424 #include <math.h>
2525 #include <stdlib.h>
2626
27-void libtetrabz_dblstep2(
27+static void libtetrabz_dblstep2(
2828 int nb,
2929 double ei[4],
3030 double **ej,
--- a/src/c/libtetrabz_fermigr.c
+++ b/src/c/libtetrabz_fermigr.c
@@ -23,7 +23,7 @@
2323 #include "libtetrabz_common.h"
2424 #include <stdlib.h>
2525
26-void libtetrabz_fermigr3(
26+static void libtetrabz_fermigr3(
2727 int ne,
2828 double *e0,
2929 double de[4],
@@ -68,7 +68,7 @@ void libtetrabz_fermigr3(
6868 }
6969 }
7070
71-void libtetrabz_fermigr2(
71+static void libtetrabz_fermigr2(
7272 int nb,
7373 int ne,
7474 double *e0,
--- a/src/c/libtetrabz_polcmplx.c
+++ b/src/c/libtetrabz_polcmplx.c
@@ -25,7 +25,7 @@
2525 #include <stdlib.h>
2626 #include <stdio.h>
2727
28-void libtetrabz_polcmplx_1234(
28+static void libtetrabz_polcmplx_1234(
2929 double g1,
3030 double g2,
3131 double g3,
@@ -75,7 +75,7 @@ void libtetrabz_polcmplx_1234(
7575 }
7676
7777
78-void libtetrabz_polcmplx_1231(
78+static void libtetrabz_polcmplx_1231(
7979 double g1,
8080 double g2,
8181 double g3,
@@ -116,7 +116,7 @@ void libtetrabz_polcmplx_1231(
116116 }
117117
118118
119-void libtetrabz_polcmplx_1233(
119+static void libtetrabz_polcmplx_1233(
120120 double g1,
121121 double g2,
122122 double g3,
@@ -164,7 +164,7 @@ void libtetrabz_polcmplx_1233(
164164 w[1] = (w3 - w2) / (2.0 * (g3 - g2));
165165 }
166166
167-void libtetrabz_polcmplx_1221(
167+static void libtetrabz_polcmplx_1221(
168168 double g1,
169169 double g2,
170170 double* w
@@ -192,7 +192,7 @@ void libtetrabz_polcmplx_1221(
192192 }
193193
194194
195-void libtetrabz_polcmplx_1222(
195+static void libtetrabz_polcmplx_1222(
196196 double g1,
197197 double g2,
198198 double* w
@@ -220,7 +220,7 @@ void libtetrabz_polcmplx_1222(
220220 }
221221
222222
223-void libtetrabz_polcmplx_1211(
223+static void libtetrabz_polcmplx_1211(
224224 double g1,
225225 double g2,
226226 double* w
@@ -247,7 +247,7 @@ void libtetrabz_polcmplx_1211(
247247 w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
248248 }
249249
250-void libtetrabz_polcmplx3(
250+static void libtetrabz_polcmplx3(
251251 int ne,
252252 double** e0,
253253 double de[4],
@@ -425,7 +425,7 @@ void libtetrabz_polcmplx3(
425425 // Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
426426 // for 0<x<1, 0<y<1-x, 0<z<1-x-y
427427
428-void libtetrabz_polcmplx2(
428+static void libtetrabz_polcmplx2(
429429 int nb,
430430 int ne,
431431 double* *e0,
--- a/src/c/libtetrabz_polstat.c
+++ b/src/c/libtetrabz_polstat.c
@@ -29,7 +29,7 @@
2929 Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
3030 for 0<x<1, 0<y<1-x, 0<z<1-x-y
3131 */
32-double libtetrabz_polstat_1234(
32+static double libtetrabz_polstat_1234(
3333 double g1,
3434 double g2,
3535 double g3,
@@ -53,7 +53,7 @@ double libtetrabz_polstat_1234(
5353 }
5454
5555
56-double libtetrabz_polstat_1231(
56+static double libtetrabz_polstat_1231(
5757 double g1,
5858 double g2,
5959 double g3,
@@ -73,7 +73,7 @@ double libtetrabz_polstat_1231(
7373 return w;
7474 }
7575
76-double libtetrabz_polstat_1233(
76+static double libtetrabz_polstat_1233(
7777 double g1,
7878 double g2,
7979 double g3,
@@ -97,7 +97,7 @@ double libtetrabz_polstat_1233(
9797 return w;
9898 }
9999
100-double libtetrabz_polstat_1221(
100+static double libtetrabz_polstat_1221(
101101 double g1,
102102 double g2,
103103 double lng1,
@@ -114,7 +114,7 @@ double libtetrabz_polstat_1221(
114114 return w;
115115 }
116116
117-double libtetrabz_polstat_1222(
117+static double libtetrabz_polstat_1222(
118118 double g1,
119119 double g2,
120120 double lng1,
@@ -131,7 +131,7 @@ double libtetrabz_polstat_1222(
131131 return w;
132132 }
133133
134-double libtetrabz_polstat_1211(
134+static double libtetrabz_polstat_1211(
135135 double g1,
136136 double g2,
137137 double lng1,
@@ -148,7 +148,7 @@ double libtetrabz_polstat_1211(
148148 return w;
149149 }
150150
151-void libtetrabz_polstat3(
151+static void libtetrabz_polstat3(
152152 double de[4],
153153 double w1[4]
154154 )
@@ -299,7 +299,7 @@ void libtetrabz_polstat3(
299299 }
300300 }
301301
302-void libtetrabz_polstat2(
302+static void libtetrabz_polstat2(
303303 int nb,
304304 double *ei1,
305305 double **ej1,
--- /dev/null
+++ b/src/c/setup.py
@@ -0,0 +1,7 @@
1+from distutils.core import setup, Extension
2+
3+setup(name='libtetrabzc',
4+ version='1.0',
5+ ext_modules=[Extension('libtetrabzc', ['libtetrabz.c'])]
6+)
7+
--- a/src/c/test.c
+++ b/src/c/test.c
@@ -23,23 +23,12 @@
2323 #include <stdio.h>
2424 #include <stdlib.h>
2525 #include <math.h>
26+#include "libtetrabz.h"
2627
2728 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(
29+
30+
31+static void test_occ(
4332 int ng[3],
4433 int nk,
4534 int nb,
@@ -80,7 +69,7 @@ void test_occ(
8069 free(eig2);
8170 }
8271
83-void test_fermieng(
72+static void test_fermieng(
8473 int ng[3],
8574 int nk,
8675 int nb,
@@ -101,6 +90,7 @@ void test_fermieng(
10190 }
10291
10392 fermieng(ng, nk, nb, bvec, eig, wght, nelec, &iteration, &ef);
93+ printf("debug %d\n", iteration);
10494
10595 for (ib = 0; ib < nb; ib++) val[ib] = 0.0;
10696 for (ik = 0; ik < nk; ik++)
@@ -118,7 +108,7 @@ void test_fermieng(
118108 }
119109
120110
121-void test_dos(
111+static void test_dos(
122112 int ng[3],
123113 int nk,
124114 int nb,
@@ -172,7 +162,7 @@ void test_dos(
172162 free(wght);
173163 }
174164
175-void test_intdos(
165+static void test_intdos(
176166 int ng[3],
177167 int nk,
178168 int nb,
@@ -227,7 +217,7 @@ void test_intdos(
227217 }
228218
229219
230-void test_dblstep(
220+static void test_dblstep(
231221 int ng[3],
232222 int nk,
233223 int nb,
@@ -288,7 +278,7 @@ void test_dblstep(
288278 }
289279
290280
291-void test_dbldelta(
281+static void test_dbldelta(
292282 int ng[3],
293283 int nk,
294284 int nb,
@@ -348,7 +338,7 @@ void test_dbldelta(
348338 }
349339
350340
351-void test_polstat(
341+static void test_polstat(
352342 int ng[3],
353343 int nk,
354344 int nb,
@@ -411,7 +401,7 @@ void test_polstat(
411401 }
412402
413403
414-void test_fermigr(
404+static void test_fermigr(
415405 int ng[3],
416406 int nk,
417407 int nb,
@@ -493,7 +483,7 @@ void test_fermigr(
493483 free(eig22);
494484 }
495485
496-void test_polcmplx(
486+static void test_polcmplx(
497487 int ng[3],
498488 int nk,
499489 int nb,
--- /dev/null
+++ b/src/c/test.py
@@ -0,0 +1,479 @@
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+import math
24+import libtetrabzc
25+
26+
27+def test_occ(ng, nk, nb, bvec, vbz, eig, mat):
28+ """
29+
30+ :return:
31+ """
32+ #
33+ eig0 = eig.copy()
34+ for ikb in range(nk * nb):
35+ eig0[ikb] += - 0.5
36+ #
37+ wght = libtetrabzc.occ(ng[0], ng[1], ng[2], nk, nb,
38+ bvec[0][0], bvec[0][1], bvec[0][2],
39+ bvec[1][0], bvec[1][1], bvec[1][2],
40+ bvec[2][0], bvec[2][1], bvec[2][2], eig0)
41+ #
42+ val = [0.0, 0.0]
43+ for ik in range(nk):
44+ for ib in range(nb):
45+ val[ib] += wght[ik * nb + ib] * mat[ik]
46+ print("# libtetrabz_occ")
47+ print(' %15.5e %15.5e' % (4.0 * math.pi / 5.0, val[0] * vbz))
48+ print(' %15.5e %15.5e' % (math.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))
49+ print("")
50+
51+
52+def test_fermieng(ng, nk, nb, bvec, vbz, eig, mat):
53+ """
54+
55+ :return:
56+ """
57+ #
58+ nelec = (4.0 * math.pi / 3.0 + math.sqrt(2.0) * math.pi / 3.0) / vbz
59+
60+ wght = libtetrabzc.fermieng(ng[0], ng[1], ng[2], nk, nb,
61+ bvec[0][0], bvec[0][1], bvec[0][2],
62+ bvec[1][0], bvec[1][1], bvec[1][2],
63+ bvec[2][0], bvec[2][1], bvec[2][2], eig, nelec)
64+ ef = wght[nk * nb]
65+ # iterations = wght[nk*nb + 1]
66+ #
67+ val = [0.0, 0.0]
68+ for ik in range(nk):
69+ for ib in range(nb):
70+ val[ib] += wght[ik * nb + ib] * mat[ik]
71+ #
72+ print("# libtetrabz_fermieng")
73+ print(" %15.5e %15.5e" % (0.5, ef))
74+ print(" %15.5e %15.5e" % (4.0 * math.pi / 5.0, val[0] * vbz))
75+ print(" %15.5e %15.5e" % (math.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))
76+ print("")
77+
78+
79+def test_dos(ng, nk, nb, bvec, vbz, eig, mat):
80+ """
81+
82+ :return:
83+ """
84+ ne = 5
85+ #
86+ e0 = [0.0, 0.0, 0.0, 0.0, 0.0]
87+ val0 = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]
88+ for ie in range(ne):
89+ e0[ie] = 0.2 * (ie + 1)
90+ val0[0][ie] = 4.0 * math.pi * e0[ie] ** 3
91+ if e0[ie] > 1.0 / math.sqrt(2.0):
92+ val0[1][ie] = math.sqrt(2.0) * math.pi * math.sqrt(-1.0 + 2.0 * e0[ie] ** 2) ** 3
93+ else:
94+ val0[1][ie] = 0.0
95+ for ie in range(ne):
96+ e0[ie] = e0[ie] ** 2 * 0.5
97+ #
98+ wght = libtetrabzc.dos(ng[0], ng[1], ng[2], nk, nb, ne,
99+ bvec[0][0], bvec[0][1], bvec[0][2],
100+ bvec[1][0], bvec[1][1], bvec[1][2],
101+ bvec[2][0], bvec[2][1], bvec[2][2], eig, e0)
102+ #
103+ val = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]
104+ for ik in range(nk):
105+ for ib in range(nb):
106+ for ie in range(ne):
107+ val[ib][ie] += wght[ik * nb * ne + ib * ne + ie] * mat[ik]
108+ #
109+ print("# libtetrabz_dos")
110+ for ib in range(nb):
111+ for ie in range(ne):
112+ print(" %15.5e %15.5e" % (val0[ib][ie], val[ib][ie] * vbz))
113+ print("")
114+
115+
116+def test_intdos(ng, nk, nb, bvec, vbz, eig, mat):
117+ """
118+
119+ :return:
120+ """
121+ ne = 5
122+ e0 = [0.0, 0.0, 0.0, 0.0, 0.0]
123+ val0 = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]
124+ #
125+ for ie in range(ne):
126+ e0[ie] = 0.2 * (ie + 1)
127+ val0[0][ie] = 4.0 * math.pi * e0[ie] ** 5 / 5.0
128+ if e0[ie] > 1.0 / math.sqrt(2.0):
129+ val0[1][ie] = math.pi * math.sqrt(-1.0 + 2.0 * e0[ie] ** 2) ** 5 / (5.0 * math.sqrt(2.0))
130+ else:
131+ val0[1][ie] = 0.0
132+ for ie in range(ne):
133+ e0[ie] = e0[ie] ** 2 * 0.5
134+ #
135+ wght = libtetrabzc.intdos(ng[0], ng[1], ng[2], nk, nb, ne,
136+ bvec[0][0], bvec[0][1], bvec[0][2],
137+ bvec[1][0], bvec[1][1], bvec[1][2],
138+ bvec[2][0], bvec[2][1], bvec[2][2], eig, e0)
139+ #
140+ val = [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]]
141+ for ik in range(nk):
142+ for ib in range(nb):
143+ for ie in range(ne):
144+ val[ib][ie] += wght[ik * nb * ne + ib * ne + ie] * mat[ik]
145+ #
146+ print("# libtetrabz_intdos")
147+ for ib in range(nb):
148+ for ie in range(ne):
149+ print(" %15.5e %15.5e" % (val0[ib][ie], val[ib][ie] * vbz))
150+ print("")
151+
152+
153+def test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat):
154+ """
155+
156+ :param bvec:
157+ :param vbz:
158+ :param eig1:
159+ :param eig2:
160+ :param mat:
161+ :return:
162+ """
163+ #
164+ eig10 = eig1.copy()
165+ eig20 = eig2.copy()
166+ for ikb in range(nk * nb):
167+ eig10[ikb] += - 0.5
168+ eig20[ikb] += - 0.5
169+
170+ wght = libtetrabzc.dblstep(ng[0], ng[1], ng[2], nk, nb,
171+ bvec[0][0], bvec[0][1], bvec[0][2],
172+ bvec[1][0], bvec[1][1], bvec[1][2],
173+ bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)
174+ #
175+ val = [[0.0, 0.0], [0.0, 0.0]]
176+ for ik in range(nk):
177+ for ib in range(nb):
178+ for jb in range(nb):
179+ val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]
180+ #
181+ print("# libtetrabz_dblstep")
182+ print(" %15.5e %15.5e" % (49.0 * math.pi / 320.0, val[0][0] * vbz))
183+ print(" %15.5e %15.5e" % (0.0, val[0][1] * vbz))
184+ print(" %15.5e %15.5e" % (math.pi * (512.0 * math.sqrt(2.0) - 319.0) / 10240.0, val[1][0] * vbz))
185+ print(" %15.5e %15.5e" % (0.0, val[1][1] * vbz))
186+ print("")
187+
188+
189+def test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat):
190+ """
191+
192+ :param bvec:
193+ :param vbz:
194+ :param eig1:
195+ :param eig2:
196+ :param mat:
197+ :return:
198+ """
199+ #
200+ eig10 = eig1.copy()
201+ eig20 = eig2.copy()
202+ for ikb in range(nk * nb):
203+ eig10[ikb] += - 0.5
204+ eig20[ikb] += - 0.5
205+
206+ wght = libtetrabzc.dbldelta(ng[0], ng[1], ng[2], nk, nb,
207+ bvec[0][0], bvec[0][1], bvec[0][2],
208+ bvec[1][0], bvec[1][1], bvec[1][2],
209+ bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)
210+ #
211+ val = [[0.0, 0.0], [0.0, 0.0]]
212+ for ik in range(nk):
213+ for ib in range(nb):
214+ for jb in range(nb):
215+ val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]
216+ #
217+ print("# libtetrabz_dbldelta")
218+ print(" %15.5e %15.5e" % (2.0 * math.pi, val[0][0] * vbz))
219+ print(" %15.5e %15.5e" % (0.0, val[0][1] * vbz))
220+ print(" %15.5e %15.5e" % (math.pi, val[1][0] * vbz))
221+ print(" %15.5e %15.5e" % (0.0, val[1][1] * vbz))
222+ print("")
223+
224+
225+def test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat):
226+ """
227+
228+ :param bvec:
229+ :param vbz:
230+ :param eig1:
231+ :param eig2:
232+ :param mat:
233+ :return:
234+ """
235+ eig10 = eig1.copy()
236+ eig20 = eig2.copy()
237+ for ikb in range(nk * nb):
238+ eig10[ikb] += - 0.5
239+ eig20[ikb] += - 0.5
240+
241+ wght = libtetrabzc.polstat(ng[0], ng[1], ng[2], nk, nb,
242+ bvec[0][0], bvec[0][1], bvec[0][2],
243+ bvec[1][0], bvec[1][1], bvec[1][2],
244+ bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20)
245+ #
246+ val = [[0.0, 0.0], [0.0, 0.0]]
247+ for ik in range(nk):
248+ for ib in range(nb):
249+ for jb in range(nb):
250+ val[ib][jb] += wght[ik * nb * nb + ib * nb + jb] * mat[ik]
251+ #
252+ print("# libtetrabz_polstat")
253+ print(" %15.5e %15.5e" % (math.pi * (68.0 + 45.0 * math.log(3.0)) / 96.0, val[0][0] * vbz))
254+ print(" %15.5e %15.5e" % (math.pi * 8.0 / 5.0, val[0][1] * vbz))
255+ print(" %15.5e %15.5e" % (math.pi * (228.0 + 22.0 * math.sqrt(2.0) - 96.0 * math.log(2.0)
256+ + 192.0 * math.log(4.0 + math.sqrt(2.0))
257+ - 3.0 * math.log(1.0 + 2.0 * math.sqrt(2.0))) / 1536.0,
258+ val[1][0] * vbz))
259+ print(" %15.5e %15.5e" % (math.pi * math.sqrt(8.0) / 5.0, val[1][1] * vbz))
260+ print("")
261+
262+
263+def test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat):
264+ """
265+
266+ :param bvec:
267+ :param vbz:
268+ :param eig1:
269+ :param eig2:
270+ :param mat:
271+ :return:
272+ """
273+ ne = 3
274+ #
275+ eig10 = eig1.copy()
276+ eig20 = eig2.copy()
277+ for ikb in range(nk * nb):
278+ eig10[ikb] += - 0.5
279+ eig20[ikb] += - 0.5
280+
281+ e0 = [0.0, 0.0, 0.0]
282+ 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]]]
283+ e0[0] = 1.0 / 3.0
284+ e0[1] = 2.0 / 3.0
285+ e0[2] = 1.0
286+ val0[0][0][0] = 4.0 * math.pi / 9.0
287+ val0[0][0][1] = 1295.0 * math.pi / 2592.0
288+ val0[0][0][2] = 15.0 * math.pi / 32.0
289+ val0[1][0][0] = 5183.0 * math.pi / 41472.0
290+ val0[1][0][1] = 4559.0 * math.pi / 41472.0
291+ val0[1][0][2] = 0.0
292+ for ib in range(nb):
293+ for ie in range(ne):
294+ val0[ib][1][ie] = 0.0
295+
296+ wght = libtetrabzc.fermigr(ng[0], ng[1], ng[2], nk, nb, ne,
297+ bvec[0][0], bvec[0][1], bvec[0][2],
298+ bvec[1][0], bvec[1][1], bvec[1][2],
299+ bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20, e0)
300+ #
301+ 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]]]
302+ for ik in range(nk):
303+ for ib in range(nb):
304+ for jb in range(nb):
305+ for ie in range(ne):
306+ val[ib][jb][ie] += wght[ik * nb * nb * ne + ib * nb * ne + jb * ne + ie] * mat[ik]
307+ #
308+ print("# libtetrabz_fermigr")
309+ for ib in range(nb):
310+ for jb in range(nb):
311+ for ie in range(ne):
312+ print(" %15.5e %15.5e" % (val0[ib][jb][ie], val[ib][jb][ie] * vbz))
313+ print("")
314+
315+
316+def test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat):
317+ """
318+
319+ :param bvec:
320+ :param vbz:
321+ :param eig1:
322+ :param eig2:
323+ :param mat:
324+ :return:
325+ """
326+ ne = 3
327+
328+ eig10 = eig1.copy()
329+ eig20 = eig2.copy()
330+ for ikb in range(nk * nb):
331+ eig10[ikb] += - 0.5
332+ eig20[ikb] += - 0.5
333+
334+ e0 = [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]
335+ e0[0][0] = -2.0
336+ e0[0][1] = 1.0
337+ e0[1][0] = 0.0
338+ e0[1][1] = 2.0
339+ e0[2][0] = 1.0
340+ e0[2][1] = - 0.5
341+ e01 = [e0[0][0], e0[0][1], e0[1][0], e0[1][1], e0[2][0], e0[2][1]]
342+
343+ 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]]],
344+ [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]]]
345+ val0[0][0][0][0] = -0.838243341280338
346+ val0[0][0][0][1] = - 0.734201894333234
347+ val0[0][0][1][0] = 0.270393588876530
348+ val0[0][0][1][1] = - 0.771908416949610
349+ val0[0][0][2][0] = 0.970996830573510
350+ val0[0][0][2][1] = 0.302792326476720
351+ val0[1][0][0][0] = -0.130765724778920
352+ val0[1][0][0][1] = - 0.087431218706638
353+ val0[1][0][1][0] = 0.030121954547245
354+ val0[1][0][1][1] = - 0.135354254293510
355+ val0[1][0][2][0] = 0.178882244951203
356+ val0[1][0][2][1] = 0.064232167683425
357+ for ie in range(3):
358+ val0[0][1][ie][0] = (8.0 * math.pi * (1.0 + 2.0 * e0[ie][0])) \
359+ / (5.0 * ((1.0 + 2.0 * e0[ie][0]) ** 2 + (2.0 * e0[ie][1]) ** 2))
360+ val0[0][1][ie][1] = (- 8.0 * math.pi * 2.0 * e0[ie][1]) \
361+ / (5.0 * ((1.0 + 2.0 * e0[ie][0]) ** 2 + (2.0 * e0[ie][1]) ** 2))
362+ val0[1][1][ie][0] = (math.sqrt(8.0) * math.pi * (1.0 + 4.0 * e0[ie][0])) \
363+ / (5.0 * ((1.0 + 4.0 * e0[ie][0]) ** 2 + (4.0 * e0[ie][1]) ** 2))
364+ val0[1][1][ie][1] = (- math.sqrt(8.0) * math.pi * 4.0 * e0[ie][1]) \
365+ / (5.0 * ((1.0 + 4.0 * e0[ie][0]) ** 2 + (4.0 * e0[ie][1]) ** 2))
366+
367+ wght = libtetrabzc.polcmplx(ng[0], ng[1], ng[2], nk, nb, ne,
368+ bvec[0][0], bvec[0][1], bvec[0][2],
369+ bvec[1][0], bvec[1][1], bvec[1][2],
370+ bvec[2][0], bvec[2][1], bvec[2][2], eig10, eig20, e01)
371+ #
372+ 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]]],
373+ [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]]]
374+ for ik in range(nk):
375+ for ib in range(nb):
376+ for jb in range(nb):
377+ for ie in range(ne):
378+ for ir in range(2):
379+ val[ib][jb][ie][ir] += wght[ik * nb * nb * ne * 2
380+ + ib * nb * ne * 2 + jb * ne * 2 + ie * 2 + ir] * mat[ik]
381+ #
382+ print("# libtetrabz_polcmplx")
383+ for ib in range(nb):
384+ for jb in range(nb):
385+ for ie in range(ne):
386+ print(" %15.5e %15.5e" % (val0[ib][jb][ie][0], val[ib][jb][ie][0] * vbz))
387+ print(" %15.5e %15.5e" % (val0[ib][jb][ie][1], val[ib][jb][ie][1] * vbz))
388+ print("")
389+
390+
391+def test():
392+ ng0 = 8
393+ ng = [ng0, ng0, ng0]
394+ nk = ng[0] * ng[1] * ng[2]
395+ nb = 2
396+ bvec = [[3.0, 0.0, 0.0],
397+ [0.0, 3.0, 0.0],
398+ [0.0, 0.0, 3.0]]
399+ vbz = bvec[0][0] * bvec[1][1] * bvec[2][2] # abs(numpy.linalg.det(bvec))
400+
401+ #
402+ eig1 = []
403+ eig2 = []
404+ mat = []
405+ #
406+ kvec = [0.0, 0.0, 0.0]
407+ for i0 in range(ng[0]):
408+ for i1 in range(ng[1]):
409+ for i2 in range(ng[2]):
410+ kvec[0] = i0 / ng[0]
411+ kvec[1] = i1 / ng[1]
412+ kvec[2] = i2 / ng[2]
413+ #
414+ for ii in range(3):
415+ if kvec[ii] > 0.5:
416+ kvec[ii] += -1
417+ # kvec[0:3] = kvec.dot(bvec)
418+ for ii in range(3):
419+ kvec[ii] *= bvec[ii][ii]
420+ #
421+ eig01 = 0.0
422+ for ii in range(3):
423+ eig01 += kvec[ii] ** 2
424+ eig01 *= 0.5
425+ eig1.append(eig01)
426+ eig1.append(eig01 + 0.25)
427+ #
428+ kvec[0] += 1.0
429+ eig02 = 0.0
430+ for ii in range(3):
431+ eig02 += kvec[ii] ** 2
432+ eig02 *= 0.5
433+ eig2.append(eig02)
434+ eig2.append(eig01 + 0.5)
435+ #
436+ #
437+ for i0 in range(ng[1]):
438+ for i1 in range(ng[1]):
439+ for i2 in range(ng[2]):
440+ kvec[0] = i0 / ng[0]
441+ kvec[1] = i1 / ng[1]
442+ kvec[2] = i2 / ng[2]
443+ #
444+ for ii in range(3):
445+ if kvec[ii] > 0.5:
446+ kvec[ii] += -1
447+ # kvec[0:3] = kvec.dot(bvec)
448+ for ii in range(3):
449+ kvec[ii] *= bvec[ii][ii]
450+ #
451+ mat0 = 0.0
452+ for ii in range(3):
453+ mat0 += kvec[ii] ** 2
454+ #
455+ mat.append(mat0)
456+ #
457+ print("# Ideal Result")
458+ #
459+ test_occ(ng, nk, nb, bvec, vbz, eig1, mat)
460+ #
461+ test_fermieng(ng, nk, nb, bvec, vbz, eig1, mat)
462+ #
463+ test_dos(ng, nk, nb, bvec, vbz, eig1, mat)
464+ #
465+ test_intdos(ng, nk, nb, bvec, vbz, eig1, mat)
466+ #
467+ test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat)
468+ #
469+ test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat)
470+ #
471+ test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat)
472+ #
473+ test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat)
474+ #
475+ test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat)
476+ #
477+
478+
479+test()
Show on old repository browser