• R/O
  • HTTP
  • SSH
  • HTTPS

python: 提交

libtetrabz python package


Commit MetaInfo

修訂78471299f204676342f69d8657e7e2f156f34c03 (tree)
時間2022-03-15 02:18:57
作者Mitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Change Summary

差異

--- /dev/null
+++ b/src/c/__init__.c
@@ -0,0 +1,31 @@
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+from .libtetrabz_dos import dos
24+from .libtetrabz_dos import intdos
25+from .libtetrabz_occ import occ
26+from .libtetrabz_occ import fermieng
27+from .libtetrabz_polstat import polstat
28+from .libtetrabz_fermigr import fermigr
29+from .libtetrabz_dbldelta import dbldelta
30+from .libtetrabz_dblstep import dblstep
31+from .libtetrabz_polcmplx import polcmplx
--- /dev/null
+++ b/src/c/libtetrabz_common.c
@@ -0,0 +1,581 @@
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+
24+
25+void libtetrabz_initialize(
26+ int *ng,
27+ double **bvec,
28+ double **wlsm,
29+ int **ikv
30+ )
31+ {
32+ /*
33+ define shortest diagonal line & define type of tetragonal
34+ */
35+ int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
36+ double bvec2[3][3], bvec3[4][3], *bnorm;
37+ double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
38+ {0.0, 1440.0, 0.0, 30.0},
39+ {30.0, 0.0, 1440.0, 0.0},
40+ {0.0, 30.0, 0.0, 1440.0} },
41+ wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
42+ {-28.0, -38.0, 7.0, 17.0},
43+ {17.0, -28.0, -38.0, 7.0},
44+ {7.0, 17.0, -28.0, -38.0}},
45+ wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
46+ {9.0, -56.0, 9.0, -46.0},
47+ {-46.0, 9.0, -56.0, 9.0},
48+ {9.0, -46.0, 9.0, -56.0}},
49+ wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
50+ {7.0, -38.0, -28.0, 17.0},
51+ {17.0, 7.0, -38.0, -28.0},
52+ {-28.0, 17.0, 7.0, -38.0}},
53+ wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
54+ {-18.0, -18.0, -18.0, 12.0},
55+ {12.0, -18.0, -18.0, -18.0},
56+ {-18.0, 12.0, -18.0, -18.0}};
57+
58+ for (i1 = 0; i1 < 3; i1++)
59+ for (i2 = 0; i2 < 3; i2++)
60+ bvec2[i1][i2] = bvec[i1][i2] / (double) ng[i1];
61+
62+ for (i1 = 0; i1 < 3; i1++) {
63+ bvec3[0][i1] = -bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
64+ bvec3[1][i1] = bvec2[0][i1] - bvec2[1][i1] + bvec2[2][i1];
65+ bvec3[2][i1] = bvec2[0][i1] + bvec2[1][i1] - bvec2[2][i1];
66+ bvec3[3][i1] = bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
67+ }
68+ /*
69+ length of delta bvec
70+ */
71+ for (i1 = 0; i1 < 4; i1++)
72+ bnorm[i1] = 0.0;
73+ for (i2 = 0; i2 < 3; i2++)
74+ bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
75+
76+ itype = 0;
77+ for (i1 = 1; i1 < 4; i1++)
78+ if (bnorm[i1] < bnorm[itype]) itype = i1;
79+ /*
80+ start & last
81+ */
82+ for (i0 = 0; i0 < 4; i0++) {
83+ ivvec0[i0] = 0;
84+ for (i1 = 0; i1 < 4; i1++)divvec[i0][i1] = 0;
85+ divvec[i0][i0] = 1;
86+ }
87+ ivvec0[itype] = 1;
88+ divvec[itype][itype] = -1;
89+ /*
90+ Corners of tetrahedron
91+ */
92+ it = 0;
93+ for (i0 = 0; i0 < 3; i0++) {
94+ for (i1 = 0; i1 < 3; i1++) {
95+ if (i1 == i0) continue;
96+ for (i2 = 0; i2 < 3; i2++) {
97+ if (i2 == i1 || i2 == i0) continue;
98+
99+ for (i3 = 0; i3 < 3; i3++) {
100+ ivvec[it][0][i3] = ivvec0[i3];
101+ ivvec[it][1][i3] = ivvec[it][0][i3] + divvec[i0][i3];
102+ ivvec[it][2][i3] = ivvec[it][1][i3] + divvec[i1][i3];
103+ ivvec[it][3][i3] = ivvec[it][2][i3] + divvec[i2][i3];
104+ }
105+
106+ it += 1;
107+ }
108+ }
109+ }
110+ /*
111+ Additional points
112+ */
113+ for (i1 = 0; i1 < 6; i1++) {
114+ for (i2 = 0; i2 < 3; i2++) {
115+ ivvec[i1][4][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][1][i2];
116+ ivvec[i1][5][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][2][i2];
117+ ivvec[i1][6][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][3][i2];
118+ ivvec[i1][7][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][0][i2];
119+
120+ ivvec[i1][8][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][2][i2];
121+ ivvec[i1][9][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][3][i2];
122+ ivvec[i1][10][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][0][i2];
123+ ivvec[i1][11][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][1][i2];
124+
125+ ivvec[i1][12][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][3][i2];
126+ ivvec[i1][13][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][0][i2];
127+ ivvec[i1][14][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][1][i2];
128+ ivvec[i1][15][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][2][i2];
129+
130+ ivvec[i1][16][i2] = ivvec[i1][3][i2] - ivvec[i1][0][i2] + ivvec[i1][1][i2];
131+ ivvec[i1][17][i2] = ivvec[i1][0][i2] - ivvec[i1][1][i2] + ivvec[i1][2][i2];
132+ ivvec[i1][18][i2] = ivvec[i1][1][i2] - ivvec[i1][2][i2] + ivvec[i1][3][i2];
133+ ivvec[i1][19][i2] = ivvec[i1][2][i2] - ivvec[i1][3][i2] + ivvec[i1][0][i2];
134+ }
135+ }
136+
137+ for (i1 = 0; i1 < 4; i1++) {
138+ for (i2 = 0; i2 < 4; i2++) {
139+ wlsm[i1][i2] = wlsm1[i1][i2] /= 1260.0;
140+ wlsm[i1][i2 + 4] = wlsm2[i1][i2] /= 1260.0;
141+ wlsm[i1][i2 + 8] = wlsm3[i1][i2] /= 1260.0;
142+ wlsm[i1][i2 + 12] = wlsm4[i1][i2] /= 1260.0;
143+ wlsm[i1][i2 + 16] = wlsm5[i1][i2] /= 1260.0;
144+ }
145+ }
146+ /*
147+ k-index for energy
148+ */
149+ nt = 0;
150+ for (i2 = 0; i2 < ng[2]; i2++) {
151+ for (i1 = 0; i1 < ng[1]; i1++) {
152+ for (i0 = 0; i0 < ng[0]; i0++) {
153+
154+ for (it = 0; it < 6; it++) {
155+
156+ for (ii = 0; ii < 20; ii++) {
157+ ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
158+ ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
159+ ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
160+ ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
161+ }
162+ nt += 1;
163+ }
164+ }
165+ }
166+ }
167+}
168+
169+
170+void libtetrabz_tsmall_a1(
171+ double *e,
172+ double e0,
173+ double *v,
174+ double **tsmall
175+) {
176+ /*
177+ Cut small tetrahedron A1
178+ */
179+ double a10, a20, a30;
180+ a10 = (e0 - e[0]) / (e[1] - e[0]);
181+ a20 = (e0 - e[0]) / (e[2] - e[0]);
182+ a30 = (e0 - e[0]) / (e[3] - e[0]);
183+
184+ *v = a10 * a20 * a30;
185+
186+ tsmall[0][0] = 1.0;
187+ tsmall[0][1] = 0.0;
188+ tsmall[0][2] = 0.0;
189+ tsmall[0][3] = 0.0;
190+ tsmall[1][0] = 1.0 - a10;
191+ tsmall[1][1] = a10;
192+ tsmall[1][2] = 0.0;
193+ tsmall[1][3] = 0.0;
194+ tsmall[2][0] = 1.0 - a20;
195+ tsmall[2][1] = 0.0;
196+ tsmall[2][2] = a20;
197+ tsmall[2][3] = 0.0;
198+ tsmall[3][0] = 1.0 - a30;
199+ tsmall[3][1] = 0.0;
200+ tsmall[3][2] = 0.0;
201+ tsmall[3][3] = a30;
202+}
203+
204+
205+void libtetrabz_tsmall_b1(
206+ double *e,
207+ double e0,
208+ double *v,
209+ double **tsmall
210+)
211+{
212+ /*
213+ Cut small tetrahedron B1
214+ :param e: numpy
215+ :return:
216+ */
217+ double a13, a20, a30;
218+ a13 = (e0 - e[3]) / (e[1] - e[3]);
219+ a20 = (e0 - e[0]) / (e[2] - e[0]);
220+ a30 = (e0 - e[0]) / (e[3] - e[0]);
221+
222+ *v = a20 * a30 * a13;
223+
224+ tsmall[0][0] = 1.0;
225+ tsmall[0][1] = 0.0;
226+ tsmall[0][2] = 0.0;
227+ tsmall[0][3] = 0.0;
228+ tsmall[1][0] = 1.0 - a20;
229+ tsmall[1][1] = 0.0;
230+ tsmall[1][2] = a20;
231+ tsmall[1][3] = 0.0;
232+ tsmall[2][0] = 1.0 - a30;
233+ tsmall[2][1] = 0.0;
234+ tsmall[2][2] = 0.0;
235+ tsmall[2][3] = a30;
236+ tsmall[3][0] = 0.0;
237+ tsmall[3][1] = a13;
238+ tsmall[3][2] = 0.0;
239+ tsmall[3][3] = 1.0 - a13;
240+}
241+
242+
243+void libtetrabz_tsmall_b2(
244+ double *e,
245+ double e0,
246+ double *v,
247+ double **tsmall
248+)
249+{
250+ /*
251+ Cut small tetrahedron B2
252+ :param e:
253+ :return:
254+ */
255+ double a21, a31;
256+ a21 = (e0 - e[1]) / (e[2] - e[1]);
257+ a31 = (e0 - e[1]) / (e[3] - e[1]);
258+
259+ *v = a21 * a31;
260+
261+ tsmall[0][0] = 1.0;
262+ tsmall[0][1] = 0.0;
263+ tsmall[0][2] = 0.0;
264+ tsmall[0][3] = 0.0;
265+ tsmall[1][0] = 0.0;
266+ tsmall[1][1] = 1.0;
267+ tsmall[1][2] = 0.0;
268+ tsmall[1][3] = 0.0;
269+ tsmall[2][0] = 0.0;
270+ tsmall[2][1] = 1.0 - a21;
271+ tsmall[2][2] = a21;
272+ tsmall[2][3] = 0.0;
273+ tsmall[3][0] = 0.0;
274+ tsmall[3][1] = 1.0 - a31;
275+ tsmall[3][2] = 0.0;
276+ tsmall[3][3] = a31;
277+}
278+
279+void libtetrabz_tsmall_b3(
280+ double *e,
281+ double e0,
282+ double *v,
283+ double **tsmall
284+)
285+{
286+ /*
287+ Cut small tetrahedron B3
288+ :param e:
289+ :return:
290+ */
291+ double a12, a20, a31;
292+ a12 = (e0 - e[2]) / (e[1] - e[2]);
293+ a20 = (e0 - e[0]) / (e[2] - e[0]);
294+ a31 = (e0 - e[1]) / (e[3] - e[1]);
295+
296+ *v = a12 * a20 * a31;
297+
298+ tsmall[0][0] = 1.0;
299+ tsmall[0][1] = 0.0;
300+ tsmall[0][2] = 0.0;
301+ tsmall[0][3] = 0.0;
302+ tsmall[1][0] = 1.0 - a20;
303+ tsmall[1][1] = 0.0;
304+ tsmall[1][2] = a20;
305+ tsmall[1][3] = 0.0;
306+ tsmall[2][0] = 0.0;
307+ tsmall[2][1] = a12;
308+ tsmall[2][2] = 1.0 - a12;
309+ tsmall[2][3] = 0.0;
310+ tsmall[3][0] = 0.0;
311+ tsmall[3][1] = 1.0 - a31;
312+ tsmall[3][2] = 0.0;
313+ tsmall[3][3] = a31;
314+}
315+
316+
317+void libtetrabz_tsmall_c1(
318+ double *e,
319+ double e0,
320+ double *v,
321+ double **tsmall
322+)
323+{
324+ /*
325+ Cut small tetrahedron C1
326+ :param e:
327+ :return:
328+ */
329+ double a32;
330+ a32 = (e0 - e[2]) / (e[3] - e[2]);
331+
332+ *v = a32;
333+
334+ tsmall[0][0] = 1.0;
335+ tsmall[0][1] = 0.0;
336+ tsmall[0][2] = 0.0;
337+ tsmall[0][3] = 0.0;
338+ tsmall[1][0] = 0.0;
339+ tsmall[1][1] = 1.0;
340+ tsmall[1][2] = 0.0;
341+ tsmall[1][3] = 0.0;
342+ tsmall[2][0] = 0.0;
343+ tsmall[2][1] = 0.0;
344+ tsmall[2][2] = 1.0;
345+ tsmall[2][3] = 0.0;
346+ tsmall[3][0] = 0.0;
347+ tsmall[3][1] = 0.0;
348+ tsmall[3][2] = 1.0 - a32;
349+ tsmall[3][3] = a32;
350+}
351+
352+
353+void libtetrabz_tsmall_c2(
354+ double *e,
355+ double e0,
356+ double *v,
357+ double **tsmall
358+)
359+{
360+ /*
361+ Cut small tetrahedron C2
362+ :param e:
363+ :return:
364+ */
365+ double a23, a31;
366+ a23 = (e0 - e[3]) / (e[2] - e[3]);
367+ a31 = (e0 - e[1]) / (e[3] - e[1]);
368+
369+ *v = a23 * a31;
370+
371+ tsmall[0][0] = 1.0;
372+ tsmall[0][1] = 0.0;
373+ tsmall[0][2] = 0.0;
374+ tsmall[0][3] = 0.0;
375+ tsmall[1][0] = 0.0;
376+ tsmall[1][1] = 1.0;
377+ tsmall[1][2] = 0.0;
378+ tsmall[1][3] = 0.0;
379+ tsmall[2][0] = 0.0;
380+ tsmall[2][1] = 1.0 - a31;
381+ tsmall[2][2] = 0.0;
382+ tsmall[2][3] = a31;
383+ tsmall[3][0] = 0.0;
384+ tsmall[3][1] = 0.0;
385+ tsmall[3][2] = a23;
386+ tsmall[3][3] = 1.0 - a23;
387+}
388+
389+void libtetrabz_tsmall_c3(
390+ double *e,
391+ double e0,
392+ double *v,
393+ double **tsmall
394+)
395+{
396+ /*
397+ Cut small tetrahedron C3
398+ :param e:
399+ :return:
400+ */
401+ double a23, a13, a30;
402+ a23 = (e0 - e[3]) / (e[2] - e[3]);
403+ a13 = (e0 - e[3]) / (e[1] - e[3]);
404+ a30 = (e0 - e[0]) / (e[3] - e[0]);
405+
406+ *v = a23 * a13 * a30;
407+
408+ tsmall[0][0] = 1.0;
409+ tsmall[0][1] = 0.0;
410+ tsmall[0][2] = 0.0;
411+ tsmall[0][3] = 0.0;
412+ tsmall[1][0] = 1.0 - a30;
413+ tsmall[1][1] = 0.0;
414+ tsmall[1][2] = 0.0;
415+ tsmall[1][3] = a30;
416+ tsmall[2][0] = 0.0;
417+ tsmall[2][1] = a13;
418+ tsmall[2][2] = 0.0;
419+ tsmall[2][3] = 1.0 - a13;
420+ tsmall[3][0] = 0.0;
421+ tsmall[3][1] = 0.0;
422+ tsmall[3][2] = a23;
423+ tsmall[3][3] = 1.0 - a23;
424+}
425+
426+
427+void libtetrabz_triangle_a1(
428+ double *e,
429+ double e0,
430+ double *v,
431+ double **tsmall
432+)
433+{
434+ /*
435+ Cut triangle A1
436+ :param e:
437+ :return:
438+ */
439+ double a10, a20, a30;
440+ a10 = (e0 - e[0]) / (e[1] - e[0]);
441+ a20 = (e0 - e[0]) / (e[2] - e[0]);
442+ a30 = (e0 - e[0]) / (e[3] - e[0]);
443+
444+ *v = 3.0 * a10 * a20 / (e[3] - e[0]);
445+
446+ tsmall[0][0] = 1.0 - a10;
447+ tsmall[0][1] = a10;
448+ tsmall[0][2] = 0.0;
449+ tsmall[0][3] = 0.0;
450+ tsmall[1][0] = 1.0 - a20;
451+ tsmall[1][1] = 0.0;
452+ tsmall[1][2] = a20;
453+ tsmall[1][3] = 0.0;
454+ tsmall[2][0] = 1.0 - a30;
455+ tsmall[2][1] = 0.0;
456+ tsmall[2][2] = 0.0;
457+ tsmall[2][3] = a30;
458+}
459+
460+void libtetrabz_triangle_b1(
461+ double *e,
462+ double e0,
463+ double *v,
464+ double **tsmall
465+) {
466+ /*
467+ Cut triangle B1
468+ :param e:
469+ :return:
470+ */
471+ double a30, a13, a20;
472+ a30 = (e0 - e[0]) / (e[3] - e[0]);
473+ a13 = (e0 - e[3]) / (e[1] - e[3]);
474+ a20 = (e0 - e[0]) / (e[2] - e[0]);
475+
476+ *v = 3.0 * a30 * a13 / (e[2] - e[0]);
477+
478+ tsmall[0][0] = 1.0 - a20;
479+ tsmall[0][1] = 0.0;
480+ tsmall[0][2] = a20;
481+ tsmall[0][3] = 0.0;
482+ tsmall[1][0] = 1.0 - a30;
483+ tsmall[1][1] = 0.0;
484+ tsmall[1][2] = 0.0;
485+ tsmall[1][3] = a30;
486+ tsmall[2][0] = 0.0;
487+ tsmall[2][1] = a13;
488+ tsmall[2][2] = 0.0;
489+ tsmall[2][3] = 1.0 - a13;
490+}
491+
492+
493+void libtetrabz_triangle_b2(
494+ double *e,
495+ double e0,
496+ double *v,
497+ double **tsmall
498+)
499+{
500+ /*
501+
502+ :param e:
503+ :returns:
504+ - x - first
505+ - y - second
506+ */
507+ double a12, a31, a20;
508+ a12 = (e0 - e[2]) / (e[1] - e[2]);
509+ a31 = (e0 - e[1]) / (e[3] - e[1]);
510+ a20 = (e0 - e[0]) / (e[2] - e[0]);
511+
512+ *v = 3.0 * a12 * a31 / (e[2] - e[0]);
513+
514+ tsmall[0][0] = 1.0 - a20;
515+ tsmall[0][1] = 0.0;
516+ tsmall[0][2] = a20;
517+ tsmall[0][3] = 0.0;
518+ tsmall[1][0] = 0.0;
519+ tsmall[1][1] = a12;
520+ tsmall[1][2] = 1.0 - a12;
521+ tsmall[1][3] = 0.0;
522+ tsmall[2][0] = 0.0;
523+ tsmall[2][1] = 1.0 - a31;
524+ tsmall[2][2] = 0.0;
525+ tsmall[2][3] = a31;
526+}
527+
528+void libtetrabz_triangle_c1(
529+ double *e,
530+ double e0,
531+ double *v,
532+ double **tsmall
533+)
534+{
535+/*
536+Cut triangle C1
537+*/
538+ double a03, a13, a23;
539+ a03 = (e0 - e[3]) / (e[0] - e[3]);
540+ a13 = (e0 - e[3]) / (e[1] - e[3]);
541+ a23 = (e0 - e[3]) / (e[2] - e[3]);
542+
543+ *v = 3.0 * a03 * a13 / (e[3] - e[2]);
544+
545+ tsmall[0][0] = a03;
546+ tsmall[0][1] = 0.0;
547+ tsmall[0][2] = 0.0;
548+ tsmall[0][3] = 1.0 - a03;
549+ tsmall[1][0] = 0.0;
550+ tsmall[1][1] = a13;
551+ tsmall[1][2] = 0.0;
552+ tsmall[1][3] = 1.0 - a13;
553+ tsmall[2][0] = 0.0;
554+ tsmall[2][1] = 0.0;
555+ tsmall[2][2] = a23;
556+ tsmall[2][3] = 1.0 - a23;
557+}
558+
559+void eig_sort(
560+ int n, //!< [in] the number of components
561+ double *key, //!< [in] Variables to be sorted [n].
562+ int *swap //!< [out] Order of index (sorted)
563+)
564+{
565+ int i, j, k;
566+
567+ for (i = 0; i < n; ++i) swap[i] = i;
568+
569+ for (i = 0; i < n - 1; ++i) {
570+ for (j = i + 1; j < n; ++j) {
571+ if (key[swap[j]] < key[swap[i]]) {
572+ /*
573+ Swap
574+ */
575+ k = swap[j];
576+ swap[j] = swap[i];
577+ swap[i] = k;
578+ }
579+ }/*for (j = i + 1; j < n; ++j)*/
580+ }/*for (i = 0; i < n - 1; ++i)*/
581+}/*eig_sort*/
--- /dev/null
+++ b/src/c/libtetrabz_common.h
@@ -0,0 +1,36 @@
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+
24+void libtetrabz_initialize(int *ng, double **bvec, double wlsm[4][20], int **ikv);
25+void libtetrabz_tsmall_a1(double *e, double e0, double *v, double tsmall[4][4]);
26+void libtetrabz_tsmall_b1(double *e, double e0, double *v, double tsmall[4][4]);
27+void libtetrabz_tsmall_b2(double *e, double e0, double *v, double tsmall[4][4]);
28+void libtetrabz_tsmall_b3(double *e, double e0, double *v, double tsmall[4][4]);
29+void libtetrabz_tsmall_c1(double *e, double e0, double *v, double tsmall[4][4]);
30+void libtetrabz_tsmall_c2(double *e, double e0, double *v, double tsmall[4][4]);
31+void libtetrabz_tsmall_c3(double *e, double e0, double *v, double tsmall[4][4]);
32+void libtetrabz_triangle_a1(double *e, double e0, double *v, double tsmall[3][4]);
33+void libtetrabz_triangle_b1(double *e, double e0, double *v, double tsmall[3][4]);
34+void libtetrabz_triangle_b2(double *e, double e0, double *v, double tsmall[3][4]);
35+void libtetrabz_triangle_c1(double *e, double e0, double *v, double tsmall[3][4]);
36+void eig_sort(int n, double *key, int *swap);
--- /dev/null
+++ b/src/c/libtetrabz_dbldelta.c
@@ -0,0 +1,249 @@
1+//
2+// Copyright (c) 2014 Mitsuaki Kawamura
3+//
4+// Permission is hereby granted, free of charge, to any person obtaining a
5+// copy of this software and associated documentation files (the
6+// "Software"), to deal in the Software without restriction, including
7+// without limitation the rights to use, copy, modify, merge, publish,
8+// distribute, sublicense, and/or sell copies of the Software, and to
9+// permit persons to whom the Software is furnished to do so, subject to
10+// the following conditions:
11+//
12+// The above copyright notice and this permission notice shall be included
13+// in all copies or substantial portions of the Software.
14+//
15+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+//
23+#include "libtetrabz_common.h"
24+#include <stdio.h>
25+#include <stdlib.h>
26+
27+void libtetrabz_dbldelta2(
28+ int nb,
29+ double **ej,
30+ double **w
31+) {
32+ /*
33+ 2nd step of tetrahedron method.
34+ */
35+ int ii, ib, indx[4];
36+ double a10, a20, a02, a12, v, e[3];
37+
38+ for (ib = 0; ib < nb; ib++) {
39+
40+ for (ii = 0; ii < 3; ii++)
41+ w[ib][ii] = 0.0;
42+
43+ for (ii = 0; ii < 3; ii++) e[ii] = ej[ii][ib];
44+ eig_sort(3, e, indx);
45+
46+ if (e[2] < 1.0e-10) {
47+ printf("Nesting ##\n");
48+ }
49+
50+ if ((e[0] < 0.0 && 0.0 <= e[1]) || (e[0] <= 0.0 && 0.0 < e[1])) {
51+
52+ a10 = (0.0 - e[0]) / (e[1] - e[0]);
53+ a20 = (0.0 - e[0]) / (e[2] - e[0]);
54+
55+ v = a10 / (e[2] - e[0]);
56+
57+ w[ib][indx[0]] = v * (2.0 - a10 - a20);
58+ w[ib][indx[1]] = v * a10;
59+ w[ib][indx[2]] = v * a20;
60+ }
61+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
62+
63+ a02 = (0.0 - e[2]) / (e[0] - e[2]);
64+ a12 = (0.0 - e[2]) / (e[1] - e[2]);
65+
66+ v = a12 / (e[2] - e[0]);
67+
68+ w[ib][indx[0]] = v * a02;
69+ w[ib][indx[1]] = v * a12;
70+ w[ib][indx[2]] = v * (2.0 - a02 - a12);
71+ }
72+ }
73+}
74+
75+void dbldelta(
76+ int* ng,
77+ int nk,
78+ int nb,
79+ double** bvec,
80+ double** eig1,
81+ double** eig2,
82+ double*** wght
83+)
84+{
85+ /*
86+ Main SUBROUTINE for Delta(E1) * Delta(E2)
87+ */
88+ int it, ik, ib, ii, jj, jb, ** ikv, indx[4];
89+ double wlsm[4][20], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[3][4], thr;
90+
91+ thr = 1.0e-10;
92+
93+ ikv = (int**)calloc(6 * nk, sizeof(int*));
94+ ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
95+ for (ik = 0; ik < 6 * nk; ik++) {
96+ ikv[ik] = ikv[0] + ik * 20;
97+ }
98+
99+ ei1 = (double**)calloc(4, sizeof(double*));
100+ ej1 = (double**)calloc(4, sizeof(double*));
101+ ei1[0] = (double*)calloc(4 * nb, sizeof(double));
102+ ej1[0] = (double*)calloc(4 * nb, sizeof(double));
103+ for (ii = 0; ii < 4; ii++) {
104+ ei1[ii] = ei1[0] + ii * nb;
105+ ej1[ii] = ej1[0] + ii * nb;
106+ }
107+
108+ w1 = (double***)malloc(nb, sizeof(double**));
109+ w1[0] = (double**)malloc(nb * nb, sizeof(double*));
110+ w1[0][0] = (double*)malloc(nb * nb * 4, sizeof(double));
111+ for (ib = 0; ib < nb; ib++) {
112+ w1[ii] = w1[0] + ib * nb;
113+ for (jb = 0; jb < nb; jb++) {
114+ w1[ib][jb] = w1[0][0] + ib * nb * 4 + jb * 4;
115+ }
116+ }
117+
118+ ej2 = (double**)calloc(3, sizeof(double*));
119+ ej2[0] = (double*)calloc(3 * nb, sizeof(double));
120+ for (ii = 0; ii < 3; ii++) {
121+ ej2[ii] = ej2[0] + ii * nb;
122+ }
123+
124+ w2 = (double**)calloc(nb, sizeof(double*));
125+ w2[0] = (double*)calloc(nb*3, sizeof(double));
126+ for (ib = 0; ib < nb; ib++) {
127+ w2[ib] = w2[0] + ib * 3;
128+ }
129+
130+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
131+
132+ for (ik = 0; ik < nk; ik++)
133+ for (ib = 0; ib < nb; ib++)
134+ for (jb = 0; jb < nb; jb++)
135+ wght[ik][ib][jb] = 0.0;
136+
137+ for (it = 0; it < 6 * nk; it++) {
138+
139+ for (ii = 0; ii < 4; ii++) {
140+ for (ib = 0; ib < nb; ib++) {
141+ ei1[ii][ib] = 0.0;
142+ ej1[ii][ib] = 0.0;
143+ for (jj = 0; jj < 20; jj++) {
144+ ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
145+ ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
146+ }
147+ }
148+ }
149+
150+ for (ib = 0; ib < nb; ib++)
151+ for (jb = 0; jb < nb; jb++)
152+ for (ii = 0; ii < 4; ii++)
153+ w1[ib][jb][ii] = 0.0;
154+
155+ for (ib = 0; ib < nb; ib++) {
156+
157+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
158+ eig_sort(4, e, indx);
159+
160+ if (e[0] < 0.0 && e[0] <= e[1]) {
161+
162+ libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
163+
164+ if (v > thr) {
165+ for (ii = 0; ii < 3; ii++)
166+ for (jj = 0; jj < 4; jj++)
167+ for (jb = 0; jb < nb; jb++)
168+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
169+ libtetrabz_dbldelta2(nb, ej2, w2);
170+ for (ii = 0; ii < 4; ii++)
171+ for (jj = 0; jj < 3; jj++)
172+ for (jb = 0; jb < nb; jb++)
173+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
174+ }
175+ }
176+ else if (e[1] < 0.0 && 0.0 <= e[2]) {
177+
178+ libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
179+
180+ if (v > thr) {
181+ for (ii = 0; ii < 3; ii++)
182+ for (jj = 0; jj < 4; jj++)
183+ for (jb = 0; jb < nb; jb++)
184+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
185+ libtetrabz_dbldelta2(nb, ej2, w2);
186+ for (ii = 0; ii < 4; ii++)
187+ for (jj = 0; jj < 3; jj++)
188+ for (jb = 0; jb < nb; jb++)
189+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
190+ }
191+ libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
192+
193+ if (v > thr) {
194+ for (ii = 0; ii < 3; ii++)
195+ for (jj = 0; jj < 4; jj++)
196+ for (jb = 0; jb < nb; jb++)
197+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
198+ libtetrabz_dbldelta2(nb, ej2, w2);
199+ for (ii = 0; ii < 4; ii++)
200+ for (jj = 0; jj < 3; jj++)
201+ for (jb = 0; jb < nb; jb++)
202+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
203+ }
204+ }
205+ else if (e[2] < 0.0 && 0.0 < e[3]) {
206+
207+ libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
208+
209+ if (v > thr) {
210+ for (ii = 0; ii < 3; ii++)
211+ for (jj = 0; jj < 4; jj++)
212+ for (jb = 0; jb < nb; jb++)
213+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
214+ libtetrabz_dbldelta2(nb, ej2, w2);
215+ for (ii = 0; ii < 4; ii++)
216+ for (jj = 0; jj < 3; jj++)
217+ for (jb = 0; jb < nb; jb++)
218+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
219+ }
220+ }
221+ else {
222+ continue;
223+ }
224+ }
225+ for (ii = 0; ii < 20; ii++)
226+ for (jj = 0; jj < 4; jj++)
227+ for (ib = 0; ib < nb; ib++)
228+ for (jb = 0; jb < nb; jb++)
229+ wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
230+ }
231+ for (ik = 0; ik < nk; ik++)
232+ for (ib = 0; ib < nb; ib++)
233+ for (jb = 0; jb < nb; jb++)
234+ wght[ik][ib][jb] /= (6.0 * (double)nk);
235+
236+ free(ikv[0]);
237+ free(ikv);
238+ free(ei1[0]);
239+ free(ei1);
240+ free(ej1[0]);
241+ free(ej1);
242+ free(ej2[0]);
243+ free(ej2);
244+ free(w1[0][0]);
245+ free(w1[0]);
246+ free(w1);
247+ free(w2[0]);
248+ free(w2);
249+}
--- /dev/null
+++ b/src/c/libtetrabz_dblstep.c
@@ -0,0 +1,309 @@
1+/*
2+ Copyright (c) 2014 Mitsuaki Kawamura
3+
4+ Permission is hereby granted, free of charge, to any person obtaining a
5+ copy of this software and associated documentation files (the
6+ "Software"), to deal in the Software without restriction, including
7+ without limitation the rights to use, copy, modify, merge, publish,
8+ distribute, sublicense, and/or sell copies of the Software, and to
9+ permit persons to whom the Software is furnished to do so, subject to
10+ the following conditions:
11+
12+ The above copyright notice and this permission notice shall be included
13+ in all copies or substantial portions of the Software.
14+
15+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
23+#include "libtetrabz_common.h"
24+#include <math.h>
25+
26+void libtetrabz_dblstep2(
27+ int nb,
28+ double *ei,
29+ double **ej,
30+ double **w1
31+) {
32+ /*
33+ Tetrahedron method for theta( - de)
34+ */
35+ int ii, jj, ib, indx[4];
36+ double v, thr, e[4], tsmall[4][4];
37+
38+ thr = 1.0e-8;
39+
40+ for (ib = 0; ib < nb; ib++) {
41+
42+ for (ii = 0; ii < 3; ii++)
43+ w1[ib][ii] = 0.0;
44+
45+ for (ii = 0; ii < 3; ii++) e[ii] = - ei[ii] + ej[ii][ib];
46+ eig_sort(4, e, indx);
47+
48+ if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
49+ /*
50+ Theta(0) = 0.5
51+ */
52+ v = 0.5;
53+ for (ii = 0; ii < 4; ii++)
54+ w1[ib][ii] = v * 0.25;
55+ }
56+ else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
57+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
58+ for (ii = 0; ii < 4; ii++)
59+ for (jj = 0; jj < 4; jj++)
60+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
61+ }
62+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
63+
64+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
65+ for (ii = 0; ii < 4; ii++)
66+ for (jj = 0; jj < 4; jj++)
67+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
68+
69+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
70+ for (ii = 0; ii < 4; ii++)
71+ for (jj = 0; jj < 4; jj++)
72+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
73+
74+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
75+ for (ii = 0; ii < 4; ii++)
76+ for (jj = 0; jj < 4; jj++)
77+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
78+ }
79+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
80+
81+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
82+ for (ii = 0; ii < 4; ii++)
83+ for (jj = 0; jj < 4; jj++)
84+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
85+
86+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
87+ for (ii = 0; ii < 4; ii++)
88+ for (jj = 0; jj < 4; jj++)
89+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
90+
91+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
92+ for (ii = 0; ii < 4; ii++)
93+ for (jj = 0; jj < 4; jj++)
94+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
95+ }
96+ else if (e[3] <= 0.0) {
97+ for (ii = 0; ii < 4; ii++)
98+ w1[ib][ii] = 0.25;
99+ }
100+ }
101+}
102+
103+void dblstep(
104+ int *ng,
105+ int nk,
106+ int nb,
107+ double **bvec,
108+ double **eig1,
109+ double **eig2,
110+ double ***wght
111+)
112+{
113+ /*
114+ Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
115+ */
116+ int it, ik, ib, ii, jj, jb, **ikv, indx[4];
117+ double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
118+
119+ ikv = (int**)calloc(6 * nk, sizeof(int*));
120+ ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
121+ for (ik = 0; ik < 6 * nk; ik++) {
122+ ikv[ik] = ikv[0] + ik * 20;
123+ }
124+
125+ ei1 = (double**)calloc(4, sizeof(double*));
126+ ej1 = (double**)calloc(4, sizeof(double*));
127+ ei1[0] = (double*)calloc(4 * nb, sizeof(double));
128+ ej1[0] = (double*)calloc(4 * nb, sizeof(double));
129+ for (ii = 0; ii < 4; ii++) {
130+ ei1[ii] = ei1[0] + ii * nb;
131+ ej1[ii] = ej1[0] + ii * nb;
132+ }
133+
134+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
135+
136+ for (ik = 0; ik < nk; ik++)
137+ for (ib = 0; ib < nb; ib++)
138+ for (jb = 0; jb < nb; jb++)
139+ wght[ik][ib][jb] = 0.0;
140+
141+ thr = 1.0e-8;
142+
143+ for(it = 0; it < 6*nk; it++) {
144+
145+ for (ii = 0; ii < 4; ii++) {
146+ for (ib = 0; ib < nb; ib++) {
147+ ei1[ii][ib] = 0.0;
148+ ej1[ii][ib] = 0.0;
149+ for (jj = 0; jj < 20; jj++) {
150+ ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
151+ ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
152+ }
153+ }
154+ }
155+
156+
157+ for (ib = 0; ib < nb; ib++) {
158+
159+ for (jb = 0; jb < nb; jb++)
160+ for (ii = 0; ii < 4; ii++)
161+ w1[ib][jb][ii] = 0.0;
162+
163+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
164+ eig_sort(4, e, indx);
165+
166+ if (e[0] <= 0.0 && 0.0 < e[1]) {
167+
168+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
169+
170+ if (v > thr) {
171+ for (ii = 0; ii < 4; ii++)
172+ for (jj = 0; jj < 4; jj++) {
173+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
174+ for (jb = 0; jb < nb; jb++)
175+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
176+ }
177+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
178+ for (ii = 0; ii < 4; ii++)
179+ for (jj = 0; jj < 4; jj++)
180+ for (jb = 0; jb < nb; jb++)
181+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
182+ }
183+ }
184+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
185+
186+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
187+
188+ if (v > thr) {
189+ for (ii = 0; ii < 4; ii++)
190+ for (jj = 0; jj < 4; jj++) {
191+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
192+ for (jb = 0; jb < nb; jb++)
193+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
194+ }
195+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
196+ for (ii = 0; ii < 4; ii++)
197+ for (jj = 0; jj < 4; jj++)
198+ for (jb = 0; jb < nb; jb++)
199+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
200+ }
201+
202+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
203+
204+ if (v > thr) {
205+ for (ii = 0; ii < 4; ii++)
206+ for (jj = 0; jj < 4; jj++) {
207+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
208+ for (jb = 0; jb < nb; jb++)
209+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
210+ }
211+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
212+ for (ii = 0; ii < 4; ii++)
213+ for (jj = 0; jj < 4; jj++)
214+ for (jb = 0; jb < nb; jb++)
215+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
216+ }
217+
218+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
219+
220+ if (v > thr) {
221+ for (ii = 0; ii < 4; ii++)
222+ for (jj = 0; jj < 4; jj++) {
223+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
224+ for (jb = 0; jb < nb; jb++)
225+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
226+ }
227+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
228+ for (ii = 0; ii < 4; ii++)
229+ for (jj = 0; jj < 4; jj++)
230+ for (jb = 0; jb < nb; jb++)
231+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
232+ }
233+ }
234+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
235+
236+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
237+
238+ if (v > thr) {
239+ for (ii = 0; ii < 4; ii++)
240+ for (jj = 0; jj < 4; jj++) {
241+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
242+ for (jb = 0; jb < nb; jb++)
243+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
244+ }
245+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
246+ for (ii = 0; ii < 4; ii++)
247+ for (jj = 0; jj < 4; jj++)
248+ for (jb = 0; jb < nb; jb++)
249+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
250+ }
251+
252+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
253+
254+ if (v > thr) {
255+ for (ii = 0; ii < 4; ii++)
256+ for (jj = 0; jj < 4; jj++) {
257+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
258+ for (jb = 0; jb < nb; jb++)
259+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
260+ }
261+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
262+ for (ii = 0; ii < 4; ii++)
263+ for (jj = 0; jj < 4; jj++)
264+ for (jb = 0; jb < nb; jb++)
265+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
266+ }
267+
268+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
269+
270+ if (v > thr) {
271+ for (ii = 0; ii < 4; ii++)
272+ for (jj = 0; jj < 4; jj++) {
273+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
274+ for (jb = 0; jb < nb; jb++)
275+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
276+ }
277+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
278+ for (ii = 0; ii < 4; ii++)
279+ for (jj = 0; jj < 4; jj++)
280+ for (jb = 0; jb < nb; jb++)
281+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
282+ }
283+ }
284+ else if (e[3] <= 0.0) {
285+ for (ii = 0; ii < 4; ii++) {
286+ ei2[ii] = ej1[ii][ib];
287+ for (jb = 0; jb < nb; jb++)
288+ ej2[ii][jb] = ej1[ii][jb];
289+ }
290+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
291+ for (ii = 0; ii < 4; ii++)
292+ for (jb = 0; jb < nb; jb++)
293+ w1[ib][jb][ii] += v * w2[jb][ii];
294+ }
295+ else {
296+ continue;
297+ }
298+ }
299+ for (ii = 0; ii < 20; ii++)
300+ for (jj = 0; jj < 4; jj++)
301+ for (ib = 0; ib < nb; ib++)
302+ for (jb = 0; jb < nb; jb++)
303+ wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
304+ }
305+ for (ik = 0; ik < nk; ik++)
306+ for (ib = 0; ib < nb; ib++)
307+ for (jb = 0; jb < nb; jb++)
308+ wght[ik][ib][jb] /= (6.0 * (double) nk);
309+}
--- /dev/null
+++ b/src/c/libtetrabz_dos.c
@@ -0,0 +1,283 @@
1+/*
2+ Copyright (c) 2022 Mitsuaki Kawamura
3+
4+ Permission is hereby granted, free of charge, to any person obtaining a
5+ copy of this software and associated documentation files (the
6+ "Software"), to deal in the Software without restriction, including
7+ without limitation the rights to use, copy, modify, merge, publish,
8+ distribute, sublicense, and/or sell copies of the Software, and to
9+ permit persons to whom the Software is furnished to do so, subject to
10+ the following conditions:
11+
12+ The above copyright notice and this permission notice shall be included
13+ in all copies or substantial portions of the Software.
14+
15+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
23+#include "libtetrabz_common.h"
24+#include <stdlib.h>
25+
26+void dos(
27+ int* ng,
28+ int nk,
29+ int nb,
30+ int ne,
31+ double** bvec,
32+ double** eig,
33+ double* e0,
34+ double*** wght) {
35+ /*
36+ Compute Dos : Delta(E - E1)
37+ */
38+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
39+ double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[3][4];
40+
41+ ikv = (int**)calloc(6 * nk, sizeof(int*));
42+ ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
43+ for (ik = 0; ik < 6 * nk; ik++) {
44+ ikv[ik] = ikv[0] + ik * 20;
45+ }
46+
47+ ei1 = (double**)calloc(4, sizeof(double*));
48+ ei1[0] = (double*)calloc(4 * nb, sizeof(double));
49+ for (ii = 0; ii < 4; ii++) {
50+ ei1[ii] = ei1[0] + ii * nb;
51+ }
52+
53+ w1 = (double***)calloc(nb, sizeof(double**));
54+ w1[0] = (double**)calloc(nb * ne, sizeof(double*));
55+ w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
56+ for (ib = 0; ib < nb; ib++) {
57+ w1[ii] = w1[0] + ib * ne;
58+ for (ie = 0; ie < ne; ie++) {
59+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
60+ }
61+ }
62+
63+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
64+
65+ for (ik = 0; ik < nk; ik++)
66+ for (ib = 0; ib < nb; ib++)
67+ for (ie = 0; ie < ne; ie++)
68+ wght[ik][ib][ie] = 0.0;
69+
70+ for (it = 0; it < 6 * nk; it++) {
71+
72+ for (ii = 0; ii < 4; ii++) {
73+ for (ib = 0; ib < nb; ib++) {
74+ ei1[ii][ib] = 0.0;
75+ for (jj = 0; jj < 20; jj++)
76+ ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
77+ }
78+ }
79+
80+ for (ib = 0; ib < nb; ib++)
81+ for (ie = 0; ie < ne; ie++)
82+ for (ii = 0; ii < 4; ii++)
83+ w1[ib][ie][ii] = 0.0;
84+
85+ for (ib = 0; ib < nb; ib++) {
86+
87+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
88+ eig_sort(4, e, indx);
89+
90+ for (ie = 0; ie < ne; ie++) {
91+
92+ if ((e[0] < e0[ie] && e0[ie] <= e[1]) || (e[0] <= e0[ie] && e0[ie] < e[1])) {
93+
94+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
95+ for (ii = 0; ii < 4; ii++)
96+ for (jj = 0; jj < 3; jj++)
97+ w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
98+
99+ }
100+ else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
101+
102+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
103+ for (ii = 0; ii < 4; ii++)
104+ for (jj = 0; jj < 3; jj++)
105+ w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
106+
107+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
108+ for (ii = 0; ii < 4; ii++)
109+ for (jj = 0; jj < 3; jj++)
110+ w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
111+ }
112+ else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
113+
114+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
115+ for (ii = 0; ii < 4; ii++)
116+ for (jj = 0; jj < 3; jj++)
117+ w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
118+ }
119+ else {
120+ continue;
121+ }
122+ for (ii = 0; ii < 20; ii++)
123+ for (jj = 0; jj < 4; jj++)
124+ for (ib = 0; ib < nb; ib++)
125+ for (ie = 0; ie < ne; ie++)
126+ wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
127+
128+ }
129+ }
130+ }
131+ for (ik = 0; ik < nk; ik++)
132+ for (ib = 0; ib < nb; ib++)
133+ for (ie = 0; ie < ne; ie++)
134+ wght[ik][ib][ie] /= (6.0 * (double)nk);
135+
136+ free(ikv[0]);
137+ free(ikv);
138+ free(ei1[0]);
139+ free(ei1);
140+ free(w1[0][0]);
141+ free(w1[0]);
142+ free(w1);
143+}
144+
145+void intdos(
146+ int* ng,
147+ int nk,
148+ int nb,
149+ int ne,
150+ double** bvec,
151+ double** eig,
152+ double* e0,
153+ double*** wght
154+) {
155+ /*
156+ Compute integrated Dos : theta(E - E1)
157+ */
158+
159+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
160+ double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[4][4];
161+
162+ ikv = (int**)calloc(6 * nk, sizeof(int*));
163+ ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
164+ for (ik = 0; ik < 6 * nk; ik++) {
165+ ikv[ik] = ikv[0] + ik * 20;
166+ }
167+
168+ ei1 = (double**)calloc(4, sizeof(double*));
169+ ei1[0] = (double*)calloc(4 * nb, sizeof(double));
170+ for (ii = 0; ii < 4; ii++) {
171+ ei1[ii] = ei1[0] + ii * nb;
172+ }
173+
174+ w1 = (double***)calloc(nb, sizeof(double**));
175+ w1[0] = (double**)calloc(nb * ne, sizeof(double*));
176+ w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
177+ for (ib = 0; ib < nb; ib++) {
178+ w1[ii] = w1[0] + ib * ne;
179+ for (ie = 0; ie < ne; ie++) {
180+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
181+ }
182+ }
183+
184+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
185+
186+ for (ik = 0; ik < nk; ik++)
187+ for (ib = 0; ib < nb; ib++)
188+ for (ie = 0; ie < ne; ie++)
189+ wght[ik][ib][ie] = 0.0;
190+
191+ for (it = 0; it < 6 * nk; it++) {
192+
193+ for (ii = 0; ii < 4; ii++) {
194+ for (ib = 0; ib < nb; ib++) {
195+ ei1[ii][ib] = 0.0;
196+ for (jj = 0; jj < 20; jj++)
197+ ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
198+ }
199+ }
200+
201+ for (ib = 0; ib < nb; ib++)
202+ for (ie = 0; ie < ne; ie++)
203+ for (ii = 0; ii < 4; ii++)
204+ w1[ib][ie][ii] = 0.0;
205+
206+ for (ib = 0; ib < nb; ib++) {
207+
208+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
209+ eig_sort(4, e, indx);
210+
211+ for (ie = 0; ie < ne; ie++) {
212+
213+ if ((e[0] <= e0[ie] && e0[ie] < e[1]) || (e[0] < e0[ie] && e0[ie] <= e[1])) {
214+ libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
215+ for (ii = 0; ii < 4; ii++)
216+ for (jj = 0; jj < 4; jj++)
217+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
218+ }
219+ else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
220+
221+ libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
222+ for (ii = 0; ii < 4; ii++)
223+ for (jj = 0; jj < 4; jj++)
224+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
225+
226+ libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
227+ for (ii = 0; ii < 4; ii++)
228+ for (jj = 0; jj < 4; jj++)
229+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
230+
231+ libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
232+ for (ii = 0; ii < 4; ii++)
233+ for (jj = 0; jj < 4; jj++)
234+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
235+
236+ }
237+ else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
238+
239+ libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
240+ for (ii = 0; ii < 4; ii++)
241+ for (jj = 0; jj < 4; jj++)
242+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
243+
244+ libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
245+ for (ii = 0; ii < 4; ii++)
246+ for (jj = 0; jj < 4; jj++)
247+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
248+
249+ libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
250+ for (ii = 0; ii < 4; ii++)
251+ for (jj = 0; jj < 4; jj++)
252+ w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
253+
254+ }
255+ else if (e[3] <= e0[ie]) {
256+ for (ii = 0; ii < 4; ii++)
257+ w1[ib][ie][ii] = 0.25;
258+ }
259+ else {
260+ continue;
261+ }
262+ }
263+ }
264+
265+ for (ii = 0; ii < 20; ii++)
266+ for (jj = 0; jj < 4; jj++)
267+ for (ib = 0; ib < nb; ib++)
268+ for (ie = 0; ie < ne; ie++)
269+ wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
270+ }
271+ for (ik = 0; ik < nk; ik++)
272+ for (ib = 0; ib < nb; ib++)
273+ for (ie = 0; ie < ne; ie++)
274+ wght[ik][ib][ie] /= (6.0 * (double)nk);
275+
276+ free(ikv[0]);
277+ free(ikv);
278+ free(ei1[0]);
279+ free(ei1);
280+ free(w1[0][0]);
281+ free(w1[0]);
282+ free(w1);
283+}
--- /dev/null
+++ b/src/c/libtetrabz_fermigr.c
@@ -0,0 +1,443 @@
1+//
2+// Copyright (c) 2014 Mitsuaki Kawamura
3+//
4+// Permission is hereby granted, free of charge, to any person obtaining a
5+// copy of this software and associated documentation files (the
6+// "Software"), to deal in the Software without restriction, including
7+// without limitation the rights to use, copy, modify, merge, publish,
8+// distribute, sublicense, and/or sell copies of the Software, and to
9+// permit persons to whom the Software is furnished to do so, subject to
10+// the following conditions:
11+//
12+// The above copyright notice and this permission notice shall be included
13+// in all copies or substantial portions of the Software.
14+//
15+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+//
23+#include "libtetrabz_common.h"
24+
25+void libtetrabz_fermigr3(
26+ int ne,
27+ double *e0,
28+ double de[4],
29+ double **w1
30+) {
31+ int ii, jj, ie, indx[3];
32+ double e[4], tsmall[3][4], v;
33+
34+ for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
35+ eig_sort(3, e, indx);
36+
37+ for (ie = 0; ie < ne; ie++) {
38+
39+ for (ii = 0; ii < 4; ii++) w1[ie][ii] = 0.0;
40+
41+ if (e[0] < e0[ie] && e0[ie] <= e[1]) {
42+
43+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
44+ for (ii = 0; ii < 4; ii++)
45+ for (jj = 0; jj < 3; jj++)
46+ w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
47+ }
48+ else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
49+
50+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
51+ for (ii = 0; ii < 4; ii++)
52+ for (jj = 0; jj < 3; jj++)
53+ w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
54+
55+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
56+ for (ii = 0; ii < 4; ii++)
57+ for (jj = 0; jj < 3; jj++)
58+ w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
59+ }
60+ else if (e[2] < e0[ie] && e0[ie] < e[3]) {
61+
62+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
63+ for (ii = 0; ii < 4; ii++)
64+ for (jj = 0; jj < 3; jj++)
65+ w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
66+ }
67+ }
68+}
69+
70+void libtetrabz_fermigr2(
71+ int nb,
72+ int ne,
73+ double *e0,
74+ double *ei1,
75+ double **ej1,
76+ double ***w1
77+) {
78+ int ib, ii, jj, ie, indx[4];
79+ double e[4], tsmall[4][4], v, de[4], thr, **w2;
80+
81+ for (ib = 0; ib < nb; ib++) {
82+
83+ for (ie = 0; ie < ne; ie++)
84+ for (ii = 0; ii < 4; ii++)
85+ w1[ib][ie][ii] = 0.0;
86+
87+ for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
88+ eig_sort(4, e, indx);
89+
90+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
91+
92+ libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
93+
94+ if (v > thr) {
95+ for (ii = 0; ii < 4; ii++) {
96+ de[ii] = 0.0;
97+ for (jj = 0; jj < 4; jj++)
98+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
99+ }
100+ libtetrabz_fermigr3(ne, e0, de, w2);
101+ for (ii = 0; ii < 4; ii++)
102+ for (jj = 0; jj < 4; jj++)
103+ for (ie = 0; ie < ne; ie++)
104+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
105+ }
106+ }
107+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
108+
109+ libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
110+
111+ if (v > thr) {
112+ for (ii = 0; ii < 4; ii++) {
113+ de[ii] = 0.0;
114+ for (jj = 0; jj < 4; jj++)
115+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
116+ }
117+ libtetrabz_fermigr3(ne, e0, de, w2);
118+ for (ii = 0; ii < 4; ii++)
119+ for (jj = 0; jj < 4; jj++)
120+ for (ie = 0; ie < ne; ie++)
121+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
122+ }
123+
124+ libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
125+
126+ if (v > thr) {
127+ for (ii = 0; ii < 4; ii++) {
128+ de[ii] = 0.0;
129+ for (jj = 0; jj < 4; jj++)
130+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
131+ }
132+ libtetrabz_fermigr3(ne, e0, de, w2);
133+ for (ii = 0; ii < 4; ii++)
134+ for (jj = 0; jj < 4; jj++)
135+ for (ie = 0; ie < ne; ie++)
136+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
137+ }
138+
139+ libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
140+
141+ if (v > thr) {
142+ for (ii = 0; ii < 4; ii++) {
143+ de[ii] = 0.0;
144+ for (jj = 0; jj < 4; jj++)
145+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
146+ }
147+ libtetrabz_fermigr3(ne, e0, de, w2);
148+ for (ii = 0; ii < 4; ii++)
149+ for (jj = 0; jj < 4; jj++)
150+ for (ie = 0; ie < ne; ie++)
151+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
152+ }
153+ }
154+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
155+
156+ libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
157+
158+ if (v > thr) {
159+ for (ii = 0; ii < 4; ii++) {
160+ de[ii] = 0.0;
161+ for (jj = 0; jj < 4; jj++)
162+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
163+ }
164+ libtetrabz_fermigr3(ne, e0, de, w2);
165+ for (ii = 0; ii < 4; ii++)
166+ for (jj = 0; jj < 4; jj++)
167+ for (ie = 0; ie < ne; ie++)
168+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
169+ }
170+
171+ libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
172+
173+ if (v > thr) {
174+ for (ii = 0; ii < 4; ii++) {
175+ de[ii] = 0.0;
176+ for (jj = 0; jj < 4; jj++)
177+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
178+ }
179+ libtetrabz_fermigr3(ne, e0, de, w2);
180+ for (ii = 0; ii < 4; ii++)
181+ for (jj = 0; jj < 4; jj++)
182+ for (ie = 0; ie < ne; ie++)
183+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
184+ }
185+
186+ libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
187+
188+ if (v > thr) {
189+ for (ii = 0; ii < 4; ii++) {
190+ de[ii] = 0.0;
191+ for (jj = 0; jj < 4; jj++)
192+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
193+ }
194+ libtetrabz_fermigr3(ne, e0, de, w2);
195+ for (ii = 0; ii < 4; ii++)
196+ for (jj = 0; jj < 4; jj++)
197+ for (ie = 0; ie < ne; ie++)
198+ w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
199+ }
200+ }
201+ else if (e[3] <= 0.0) {
202+ for (ii = 0; ii < 4; ii++)
203+ de[ii] = ej1[ii][ib] - ei1[ii];
204+ libtetrabz_fermigr3(ne, e0, de, w2);
205+ for (ii = 0; ii < 4; ii++)
206+ for (ie = 0; ie < ne; ie++)
207+ w1[ib][ie][ii] += w2[ie][ii];
208+ }
209+ }
210+}
211+
212+void fermigr(
213+ int *ng,
214+ int nk,
215+ int nb,
216+ int ne,
217+ double **bvec,
218+ double **eig1,
219+ double **eig2,
220+ double *e0,
221+ double ****wght
222+)
223+{
224+ /*
225+ Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
226+ */
227+ int it, ik, ib, ii, jj, jb, ie, **ikv, indx[4];
228+ double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ****w1, ***w2, v, tsmall[4][4], thr;
229+
230+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
231+
232+ for (ik = 0; ik < nk; ik++)
233+ for (ib = 0; ib < nb; ib++)
234+ for (jb = 0; jb < nb; jb++)
235+ for (ie = 0; ie < ne; ie++)
236+ wght[ik][ib][jb][ie] = 0.0;
237+
238+ thr = 1.0e-10;
239+
240+ for(it = 0; it < 6*nk; it++) {
241+
242+ for (ii = 0; ii < 4; ii++) {
243+ for (ib = 0; ib < nb; ib++) {
244+ ei1[ii][ib] = 0.0;
245+ ej1[ii][ib] = 0.0;
246+ }
247+ }
248+ for (jj = 0; jj < 20; jj++) {
249+ for (ib = 0; ib < nb; ib++) {
250+ for (ii = 0; ii < 4; ii++) {
251+ ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
252+ ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
253+ }
254+ }
255+ }
256+
257+ for (ib = 0; ib < nb; ib++) {
258+
259+ for (jb = 0; jb < nb; jb++)
260+ for (ii = 0; ii < 4; ii++)
261+ for (ie = 0; ie < ne; ie++)
262+ w1[ib][jb][ie][ii] = 0.0;
263+
264+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
265+ eig_sort(4, e, indx);
266+
267+ if (e[0] <= 0.0 && 0.0 < e[1]) {
268+
269+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
270+
271+ if (v > thr) {
272+ for (ii = 0; ii < 4; ii++) {
273+ ei2[ii] = 0.0;
274+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
275+ for (jj = 0; jj < 4; jj++) {
276+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
277+ for (jb = 0; jb < nb; jb++)
278+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
279+ }
280+ }
281+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
282+ for (ii = 0; ii < 4; ii++)
283+ for (jj = 0; jj < 4; jj++)
284+ for (jb = 0; jb < nb; jb++)
285+ for (ie = 0; ie < ne; ie++)
286+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
287+ }
288+ }
289+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
290+
291+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
292+
293+ if (v > thr) {
294+ for (ii = 0; ii < 4; ii++) {
295+ ei2[ii] = 0.0;
296+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
297+ for (jj = 0; jj < 4; jj++) {
298+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
299+ for (jb = 0; jb < nb; jb++)
300+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
301+ }
302+ }
303+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
304+ for (ii = 0; ii < 4; ii++)
305+ for (jj = 0; jj < 4; jj++)
306+ for (jb = 0; jb < nb; jb++)
307+ for (ie = 0; ie < ne; ie++)
308+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
309+ }
310+
311+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
312+
313+ if (v > thr) {
314+ for (ii = 0; ii < 4; ii++) {
315+ ei2[ii] = 0.0;
316+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
317+ for (jj = 0; jj < 4; jj++) {
318+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
319+ for (jb = 0; jb < nb; jb++)
320+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
321+ }
322+ }
323+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
324+ for (ii = 0; ii < 4; ii++)
325+ for (jj = 0; jj < 4; jj++)
326+ for (jb = 0; jb < nb; jb++)
327+ for (ie = 0; ie < ne; ie++)
328+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
329+ }
330+
331+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
332+
333+ if (v > thr) {
334+ for (ii = 0; ii < 4; ii++) {
335+ ei2[ii] = 0.0;
336+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
337+ for (jj = 0; jj < 4; jj++) {
338+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
339+ for (jb = 0; jb < nb; jb++)
340+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
341+ }
342+ }
343+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
344+ for (ii = 0; ii < 4; ii++)
345+ for (jj = 0; jj < 4; jj++)
346+ for (jb = 0; jb < nb; jb++)
347+ for (ie = 0; ie < ne; ie++)
348+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
349+ }
350+ }
351+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
352+
353+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
354+
355+ if (v > thr) {
356+ for (ii = 0; ii < 4; ii++) {
357+ ei2[ii] = 0.0;
358+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
359+ for (jj = 0; jj < 4; jj++) {
360+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
361+ for (jb = 0; jb < nb; jb++)
362+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
363+ }
364+ }
365+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
366+ for (ii = 0; ii < 4; ii++)
367+ for (jj = 0; jj < 4; jj++)
368+ for (jb = 0; jb < nb; jb++)
369+ for (ie = 0; ie < ne; ie++)
370+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
371+ }
372+
373+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
374+
375+ if (v > thr) {
376+ for (ii = 0; ii < 4; ii++) {
377+ ei2[ii] = 0.0;
378+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
379+ for (jj = 0; jj < 4; jj++) {
380+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
381+ for (jb = 0; jb < nb; jb++)
382+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
383+ }
384+ }
385+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
386+ for (ii = 0; ii < 4; ii++)
387+ for (jj = 0; jj < 4; jj++)
388+ for (jb = 0; jb < nb; jb++)
389+ for (ie = 0; ie < ne; ie++)
390+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
391+ }
392+
393+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
394+
395+ if (v > thr) {
396+ for (ii = 0; ii < 4; ii++) {
397+ ei2[ii] = 0.0;
398+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
399+ for (jj = 0; jj < 4; jj++) {
400+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
401+ for (jb = 0; jb < nb; jb++)
402+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
403+ }
404+ }
405+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
406+ for (ii = 0; ii < 4; ii++)
407+ for (jj = 0; jj < 4; jj++)
408+ for (jb = 0; jb < nb; jb++)
409+ for (ie = 0; ie < ne; ie++)
410+ w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
411+ }
412+ }
413+ else if (e[3] <= 0.0) {
414+
415+ for (ii = 0; ii < 4; ii++) {
416+ ei2[ii] = ej1[ii][ib];
417+ for (jb = 0; jb < nb; jb++)
418+ ej2[ii][jb] = ej1[ii][jb];
419+ }
420+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
421+
422+ for (ii = 0; ii < 4; ii++)
423+ for (jb = 0; jb < nb; jb++)
424+ for (ie = 0; ie < ne; ie++)
425+ w1[ib][jb][ie][ii] += w2[jb][ie][ii];
426+ }
427+ else {
428+ continue;
429+ }
430+ }
431+ for (ii = 0; ii < 20; ii++)
432+ for (jj = 0; jj < 4; jj++)
433+ for (ib = 0; ib < nb; ib++)
434+ for (jb = 0; jb < nb; jb++)
435+ for (ie = 0; ie < ne; ie++)
436+ wght[ikv[it][ii]][ib][jb][ie] += w1[ib][jb][ie][jj] * wlsm[jj][ii];
437+ }
438+ for (ik = 0; ik < nk; ik++)
439+ for (ib = 0; ib < nb; ib++)
440+ for (jb = 0; jb < nb; jb++)
441+ for (ie = 0; ie < ne; ie++)
442+ wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
443+}
--- /dev/null
+++ b/src/c/libtetrabz_occ.c
@@ -0,0 +1,207 @@
1+//
2+// Copyright (c) 2014 Mitsuaki Kawamura
3+//
4+// Permission is hereby granted, free of charge, to any person obtaining a
5+// copy of this software and associated documentation files (the
6+// "Software"), to deal in the Software without restriction, including
7+// without limitation the rights to use, copy, modify, merge, publish,
8+// distribute, sublicense, and/or sell copies of the Software, and to
9+// permit persons to whom the Software is furnished to do so, subject to
10+// the following conditions:
11+//
12+// The above copyright notice and this permission notice shall be included
13+// in all copies or substantial portions of the Software.
14+//
15+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+//
23+#include <math.h>
24+#include <stdlib.h>
25+#include "libtetrabz_common.h"
26+
27+void occ(
28+ int* ng,
29+ int nk,
30+ int nb,
31+ double** bvec,
32+ double** eig,
33+ double** wght
34+) {
35+ /*
36+ Main SUBROUTINE for occupation : Theta(EF - E1)
37+ */
38+ int it, ik, ib, ii, jj, ** ikv, indx[4];
39+ double wlsm[4][20], ** ei1, e[4], ** w1, v, tsmall[4][4];
40+
41+ ikv = (int**)calloc(6 * nk, sizeof(int*));
42+ ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
43+ for (ik = 0; ik < 6 * nk; ik++) {
44+ ikv[ik] = ikv[0] + ik * 20;
45+ }
46+
47+ ei1 = (double**)calloc(4, sizeof(double*));
48+ ei1[0] = (double*)calloc(4 * nb, sizeof(double));
49+ for (ii = 0; ii < 4; ii++) {
50+ ei1[ii] = ei1[0] + ii * nb;
51+ }
52+
53+ w1 = (double**)calloc(nb, sizeof(double*));
54+ w1[0] = (double*)calloc(nb * 4, sizeof(double));
55+ for (ib = 0; ib < nb; ib++) {
56+ w1[ii] = w1[0] + ib * 4;
57+ }
58+
59+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
60+
61+ for (ik = 0; ik < nk; ik++)
62+ for (ib = 0; ib < nb; ib++)
63+ wght[ik][ib] = 0.0;
64+
65+ for (it = 0; it < 6 * nk; it++) {
66+
67+ for (ii = 0; ii < 4; ii++) {
68+ for (ib = 0; ib < nb; ib++) {
69+ ei1[ii][ib] = 0.0;
70+ for (jj = 0; jj < 20; jj++)
71+ ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
72+ }
73+ }
74+
75+ for (ib = 0; ib < nb; ib++)
76+ for (ii = 0; ii < 4; ii++)
77+ w1[ib][ii] = 0.0;
78+
79+ for (ib = 0; ib < nb; ib++) {
80+
81+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
82+ eig_sort(4, e, indx);
83+
84+ if (e[0] <= 0.0 && 0.0 < e[1]) {
85+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
86+ for (ii = 0; ii < 4; ii++)
87+ for (jj = 0; jj < 4; jj++)
88+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
89+ }
90+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
91+
92+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
93+ for (ii = 0; ii < 4; ii++)
94+ for (jj = 0; jj < 4; jj++)
95+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
96+
97+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
98+ for (ii = 0; ii < 4; ii++)
99+ for (jj = 0; jj < 4; jj++)
100+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
101+
102+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
103+ for (ii = 0; ii < 4; ii++)
104+ for (jj = 0; jj < 4; jj++)
105+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
106+ }
107+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
108+
109+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
110+ for (ii = 0; ii < 4; ii++)
111+ for (jj = 0; jj < 4; jj++)
112+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
113+
114+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
115+ for (ii = 0; ii < 4; ii++)
116+ for (jj = 0; jj < 4; jj++)
117+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
118+
119+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
120+ for (ii = 0; ii < 4; ii++)
121+ for (jj = 0; jj < 4; jj++)
122+ w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
123+ }
124+ else if (e[3] <= 0.0) {
125+ for (ii = 0; ii < 4; ii++)
126+ w1[ib][ii] = 0.25;
127+ }
128+ else {
129+ continue;
130+ }
131+ }
132+ for (ii = 0; ii < 20; ii++)
133+ for (jj = 0; jj < 4; jj++)
134+ for (ib = 0; ib < nb; ib++)
135+ wght[ikv[it][ii]][ib] += w1[ib][jj] * wlsm[jj][ii];
136+ }
137+ for (ik = 0; ik < nk; ik++)
138+ for (ib = 0; ib < nb; ib++)
139+ wght[ik][ib] /= (6.0 * (double)nk);
140+
141+ free(ikv[0]);
142+ free(ikv);
143+ free(ei1[0]);
144+ free(ei1);
145+ free(w1[0]);
146+ free(w1);
147+}
148+
149+void fermieng(
150+ int *ng,
151+ int nk,
152+ int nb,
153+ double **bvec,
154+ double **eig,
155+ double **wght,
156+ double nelec,
157+ int *iteration,
158+ double *ef
159+) {
160+ /*
161+ Calculate Fermi energy
162+ */
163+ int maxiter, ik, ib;
164+ double eps, elw, eup, **eig2, sumkmid;
165+
166+ maxiter = 300;
167+ eps = 1.0e-10;
168+
169+ elw = eig[0][0];
170+ eup = eig[0][0];
171+ for (ik = 0; ik < nk; ik++) {
172+ for (ib = 0; ib < nb; ib++) {
173+ if (elw > eig[ik][ib]) elw = eig[ik][ib];
174+ if (eup < eig[ik][ib]) eup = eig[ik][ib];
175+ }
176+ }
177+ /*
178+ Bisection method
179+ */
180+ for (*iteration = 0; *iteration < maxiter; *iteration++) {
181+
182+ *ef = (eup + elw) * 0.5;
183+ /*
184+ Calc. # of electrons
185+ */
186+ for (ik = 0; ik < nk; ik++)
187+ for (ib = 0; ib < nb; ib++)
188+ eig2[ik][ib] = eig[ik][ib] - *ef;
189+ occ(ng, nk, nb, bvec, eig2, wght);
190+
191+ sumkmid = 0.0;
192+ for (ik = 0; ik < nk; ik++) {
193+ for (ib = 0; ib < nb; ib++) {
194+ sumkmid += wght[ik][ib];
195+ }
196+ }
197+ /*
198+ convergence check
199+ */
200+ if (fabs(sumkmid - nelec) < eps)
201+ return;
202+ else if (sumkmid < nelec)
203+ elw = *ef;
204+ else
205+ eup = *ef;
206+ }
207+}
--- /dev/null
+++ b/src/c/libtetrabz_polcmplx.c
@@ -0,0 +1,823 @@
1+/*
2+ Copyright (c) 2014 Mitsuaki Kawamura
3+
4+ Permission is hereby granted, free of charge, to any person obtaining a
5+ copy of this software and associated documentation files (the
6+ "Software"), to deal in the Software without restriction, including
7+ without limitation the rights to use, copy, modify, merge, publish,
8+ distribute, sublicense, and/or sell copies of the Software, and to
9+ permit persons to whom the Software is furnished to do so, subject to
10+ the following conditions:
11+
12+ The above copyright notice and this permission notice shall be included
13+ in all copies or substantial portions of the Software.
14+
15+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
23+#include "libtetrabz_common.h"
24+#include <math.h>
25+
26+void polcmplx(
27+ int *ng,
28+ int nk,
29+ int nb,
30+ int ne,
31+ double **bvec,
32+ double **eig1,
33+ double **eig2,
34+ double **e0,
35+ double *****wght) {
36+ /*
37+ Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
38+ */
39+ int it, ik, ie, ib, ii, jj, kk, jb, **ikv, indx[4];
40+ double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], *****w1, ****w2, v, tsmall[4][4], thr;
41+
42+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
43+
44+ for (ik = 0; ik < nk; ik++)
45+ for (ib = 0; ib < nb; ib++)
46+ for (jb = 0; jb < nb; jb++)
47+ for (ie = 0; ie < ne; ie++)
48+ for (jj = 0; jj < 2; jj++)
49+ wght[ik][ib][jb][ie][jj] = 0.0;
50+
51+ thr = 1.0e-8;
52+
53+ for (it = 0; it < 6 * nk; it++) {
54+
55+ for (ii = 0; ii < 4; ii++) {
56+ for (ib = 0; ib < nb; ib++) {
57+ ei1[ii][ib] = 0.0;
58+ ej1[ii][ib] = 0.0;
59+ }
60+ }
61+ for (jj = 0; jj < 20; jj++) {
62+ for (ib = 0; ib < nb; ib++) {
63+ for (ii = 0; ii < 4; ii++) {
64+ ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
65+ ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
66+ }
67+ }
68+ }
69+
70+ for (ib = 0; ib < nb; ib++) {
71+
72+ for (jb = 0; jb < nb; jb++)
73+ for (ii = 0; ii < 4; ii++)
74+ for (ie = 0; ie < ne; ie++)
75+ for (jj = 0; jj < 2; jj++)
76+ w1[ib][jb][ie][ii][jj] = 0.0;
77+
78+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
79+ eig_sort(4, e, indx);
80+
81+ if (e[0] <= 0.0 && 0.0 < e[1]) {
82+
83+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
84+
85+ if (v > thr) {
86+ for (ii = 0; ii < 4; ii++) {
87+ ei2[ii] = 0.0;
88+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
89+ for (jj = 0; jj < 4; jj++) {
90+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
91+ for (jb = 0; jb < nb; jb++)
92+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
93+ }
94+ }
95+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
96+ for (ii = 0; ii < 4; ii++)
97+ for (jj = 0; jj < 4; jj++)
98+ for (jb = 0; jb < nb; jb++)
99+ for (ie = 0; ie < ne; ie++)
100+ for (kk = 0; kk < 4; kk++)
101+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
102+ }
103+ }
104+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
105+
106+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
107+
108+ if (v > thr) {
109+ for (ii = 0; ii < 4; ii++) {
110+ ei2[ii] = 0.0;
111+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
112+ for (jj = 0; jj < 4; jj++) {
113+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
114+ for (jb = 0; jb < nb; jb++)
115+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
116+ }
117+ }
118+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
119+ for (ii = 0; ii < 4; ii++)
120+ for (jj = 0; jj < 4; jj++)
121+ for (jb = 0; jb < nb; jb++)
122+ for (ie = 0; ie < ne; ie++)
123+ for (kk = 0; kk < 4; kk++)
124+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
125+ }
126+
127+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
128+
129+ if (v > thr) {
130+ for (ii = 0; ii < 4; ii++) {
131+ ei2[ii] = 0.0;
132+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
133+ for (jj = 0; jj < 4; jj++) {
134+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
135+ for (jb = 0; jb < nb; jb++)
136+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
137+ }
138+ }
139+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
140+ for (ii = 0; ii < 4; ii++)
141+ for (jj = 0; jj < 4; jj++)
142+ for (jb = 0; jb < nb; jb++)
143+ for (ie = 0; ie < ne; ie++)
144+ for (kk = 0; kk < 4; kk++)
145+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
146+ }
147+
148+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
149+
150+ if (v > thr) {
151+ for (ii = 0; ii < 4; ii++) {
152+ ei2[ii] = 0.0;
153+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
154+ for (jj = 0; jj < 4; jj++) {
155+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
156+ for (jb = 0; jb < nb; jb++)
157+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
158+ }
159+ }
160+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
161+ for (ii = 0; ii < 4; ii++)
162+ for (jj = 0; jj < 4; jj++)
163+ for (jb = 0; jb < nb; jb++)
164+ for (ie = 0; ie < ne; ie++)
165+ for (kk = 0; kk < 4; kk++)
166+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
167+ }
168+ }
169+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
170+
171+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
172+
173+ if (v > thr) {
174+ for (ii = 0; ii < 4; ii++) {
175+ ei2[ii] = 0.0;
176+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
177+ for (jj = 0; jj < 4; jj++) {
178+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
179+ for (jb = 0; jb < nb; jb++)
180+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
181+ }
182+ }
183+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
184+ for (ii = 0; ii < 4; ii++)
185+ for (jj = 0; jj < 4; jj++)
186+ for (jb = 0; jb < nb; jb++)
187+ for (ie = 0; ie < ne; ie++)
188+ for (kk = 0; kk < 4; kk++)
189+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
190+ }
191+
192+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
193+
194+ if (v > thr) {
195+ for (ii = 0; ii < 4; ii++) {
196+ ei2[ii] = 0.0;
197+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
198+ for (jj = 0; jj < 4; jj++) {
199+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
200+ for (jb = 0; jb < nb; jb++)
201+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
202+ }
203+ }
204+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
205+ for (ii = 0; ii < 4; ii++)
206+ for (jj = 0; jj < 4; jj++)
207+ for (jb = 0; jb < nb; jb++)
208+ for (ie = 0; ie < ne; ie++)
209+ for (kk = 0; kk < 4; kk++)
210+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
211+ }
212+
213+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
214+
215+ if (v > thr) {
216+ for (ii = 0; ii < 4; ii++) {
217+ ei2[ii] = 0.0;
218+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
219+ for (jj = 0; jj < 4; jj++) {
220+ ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
221+ for (jb = 0; jb < nb; jb++)
222+ ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
223+ }
224+ }
225+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
226+ for (ii = 0; ii < 4; ii++)
227+ for (jj = 0; jj < 4; jj++)
228+ for (jb = 0; jb < nb; jb++)
229+ for (ie = 0; ie < ne; ie++)
230+ for (kk = 0; kk < 4; kk++)
231+ w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
232+ }
233+ }
234+ else if (e[3] <= 0.0) {
235+
236+ for (ii = 0; ii < 4; ii++) {
237+ ei2[ii] = ej1[ii][ib];
238+ for (jb = 0; jb < nb; jb++) ej2[ii][jb] = ej1[ii][jb];
239+ }
240+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
241+ for (ii = 0; ii < 4; ii++)
242+ for (jb = 0; jb < nb; jb++)
243+ for (ie = 0; ie < ne; ie++)
244+ for (kk = 0; kk < 4; kk++)
245+ w1[ib][jb][ie][kk][ii] += w2[jb][ie][kk][ii];
246+ }
247+ else {
248+ continue;
249+ }
250+ }
251+ for (ii = 0; ii < 20; ii++)
252+ for (jj = 0; jj < 4; jj++)
253+ for (ib = 0; ib < nb; ib++)
254+ for (jb = 0; jb < nb; jb++)
255+ for (ie = 0; ie < ne; ie++)
256+ for (kk = 0; kk < 2; kk++)
257+ wght[ikv[it][ii]][ib][jb][ie][kk] += w1[ib][jb][ie][kk][jj] * wlsm[jj][ii];
258+ }
259+ for (ik = 0; ik < nk; ik++)
260+ for (ib = 0; ib < nb; ib++)
261+ for (jb = 0; jb < nb; jb++)
262+ for (ie = 0; ie < ne; ie++)
263+ for (ii = 0; ii < 2; ii++)
264+ wght[ik][ib][jb][ie][ii] /= (6.0 * (double) nk);
265+}
266+
267+void libtetrabz_polcmplx2(
268+ int nb,
269+ int ne,
270+ double *e0,
271+ double *ei1,
272+ double **ej1,
273+ double ****w1
274+) {
275+ /*
276+ Tetrahedron method for theta( - E2)
277+ */
278+ int ib, ii, jj, kk, ie, indx[4];
279+ double e[4], tsmall[4][4], v, de[4], thr, ***w2;
280+
281+ thr = 1.0e-8;
282+
283+ for (ib = 0; ib < nb; ib++) {
284+
285+ for (ie = 0; ie < ne; ie++)
286+ for (ii = 0; ii < 4; ii++)
287+ for (jj = 0; jj < 2; jj++)
288+ w1[ib][ie][jj][ii] = 0.0;
289+
290+ for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
291+ eig_sort(4, e, indx);
292+
293+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
294+
295+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
296+
297+ if (v > thr) {
298+ for (ii = 0; ii < 4; ii++) {
299+ de[ii] = 0.0;
300+ for (jj = 0; jj < 4; jj++)
301+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
302+ }
303+ libtetrabz_polcmplx3(ne, e0, de, w2);
304+ for (ii = 0; ii < 4; ii++)
305+ for (jj = 0; jj < 4; jj++)
306+ for (ie = 0; ie < ne; ie++)
307+ for (kk = 0; kk < 2; kk++)
308+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
309+ }
310+ }
311+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
312+
313+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
314+
315+ if (v > thr) {
316+ for (ii = 0; ii < 4; ii++) {
317+ de[ii] = 0.0;
318+ for (jj = 0; jj < 4; jj++)
319+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
320+ }
321+ libtetrabz_polcmplx3(ne, e0, de, w2);
322+ for (ii = 0; ii < 4; ii++)
323+ for (jj = 0; jj < 4; jj++)
324+ for (ie = 0; ie < ne; ie++)
325+ for (kk = 0; kk < 2; kk++)
326+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
327+ }
328+
329+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
330+
331+ if (v > thr) {
332+ for (ii = 0; ii < 4; ii++) {
333+ de[ii] = 0.0;
334+ for (jj = 0; jj < 4; jj++)
335+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
336+ }
337+ libtetrabz_polcmplx3(ne, e0, de, w2);
338+ for (ii = 0; ii < 4; ii++)
339+ for (jj = 0; jj < 4; jj++)
340+ for (ie = 0; ie < ne; ie++)
341+ for (kk = 0; kk < 2; kk++)
342+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
343+ }
344+
345+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
346+
347+ if (v > thr) {
348+ for (ii = 0; ii < 4; ii++) {
349+ de[ii] = 0.0;
350+ for (jj = 0; jj < 4; jj++)
351+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
352+ }
353+ libtetrabz_polcmplx3(ne, e0, de, w2);
354+ for (ii = 0; ii < 4; ii++)
355+ for (jj = 0; jj < 4; jj++)
356+ for (ie = 0; ie < ne; ie++)
357+ for (kk = 0; kk < 2; kk++)
358+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
359+ }
360+ }
361+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
362+
363+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
364+
365+ if (v > thr) {
366+ for (ii = 0; ii < 4; ii++) {
367+ de[ii] = 0.0;
368+ for (jj = 0; jj < 4; jj++)
369+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
370+ }
371+ libtetrabz_polcmplx3(ne, e0, de, w2);
372+ for (ii = 0; ii < 4; ii++)
373+ for (jj = 0; jj < 4; jj++)
374+ for (ie = 0; ie < ne; ie++)
375+ for (kk = 0; kk < 2; kk++)
376+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
377+ }
378+
379+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
380+
381+ if (v > thr) {
382+ for (ii = 0; ii < 4; ii++) {
383+ de[ii] = 0.0;
384+ for (jj = 0; jj < 4; jj++)
385+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
386+ }
387+ libtetrabz_polcmplx3(ne, e0, de, w2);
388+ for (ii = 0; ii < 4; ii++)
389+ for (jj = 0; jj < 4; jj++)
390+ for (ie = 0; ie < ne; ie++)
391+ for (kk = 0; kk < 2; kk++)
392+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
393+ }
394+
395+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
396+
397+ if (v > thr) {
398+ for (ii = 0; ii < 4; ii++) {
399+ de[ii] = 0.0;
400+ for (jj = 0; jj < 4; jj++)
401+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
402+ }
403+ libtetrabz_polcmplx3(ne, e0, de, w2);
404+ for (ii = 0; ii < 4; ii++)
405+ for (jj = 0; jj < 4; jj++)
406+ for (ie = 0; ie < ne; ie++)
407+ for (kk = 0; kk < 2; kk++)
408+ w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
409+ }
410+ }
411+ else if (e[3] <= 0.0) {
412+ for (ii = 0; ii < 4; ii++) {
413+ de[ii] = tsmall[ii][jj] * (ej1[ii][ib] - ei1[ii]);
414+ }
415+ libtetrabz_polcmplx3(ne, e0, de, w2);
416+ for (ii = 0; ii < 4; ii++)
417+ for (ie = 0; ie < ne; ie++)
418+ for (kk = 0; kk < 2; kk++)
419+ w1[ib][ie][kk][ii] += w2[ie][kk][ii];
420+ }
421+ }
422+}
423+
424+
425+void libtetrabz_polcmplx3(
426+ int ne,
427+ double **e0,
428+ double de[4],
429+ double ***w1
430+) {
431+ /*
432+ Tetrahedron method for delta(om - ep + e)
433+ */
434+ int ii, jj, indx[4], ie;
435+ double e[4], x[4], thr, w2[4][2];
436+
437+ for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
438+ eig_sort(4, e, indx);
439+
440+ for (ie = 0; ie < ne; ie++) {
441+ /*
442+ I am not sure which one is better.
443+ The former is more stable.
444+ The latter is more accurate ?
445+ */
446+ for (ii = 0; ii < 4; ii++)
447+ for (jj = 0; jj < 2; jj++)
448+ w1[ie][jj][ii] = 0.25 / (de[ii] + e0[ie][jj]);
449+
450+ continue;
451+
452+ for (ii = 0; ii < 4; ii++)
453+ x[ii] = (e[ii] + e0[ie][0]) / e0[ie][1];
454+ /* thr = maxval(de(1:4)) * 1d-3*/
455+ thr = fabs(x[0]);
456+ for(ii=0;ii<4;ii++)
457+ if(thr < fabs(x[ii])) thr = fabs(x[ii]);
458+ thr = max(1.0e-3, thr * 1.0e-2);
459+
460+ if (fabs(x[3] - x[2]) < thr) {
461+ if (fabs(x[3] - x[1]) < thr) {
462+ if (fabs(x[3] - x[0]) < thr) {
463+ /*
464+ e[3] = e[2] = e[1] = e[0]
465+ */
466+ w2[3][0] = 0.25 * x[3] / (1.0 + x[3] *x[3]);
467+ w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
468+ for (ii = 0; ii < 2; ii++) {
469+ w2[2][ii] = w2[3][ii];
470+ w2[1][ii] = w2[3][ii];
471+ w2[0][ii] = w2[3][ii];
472+ }
473+ }
474+ else {
475+ /*
476+ e[3] = e[2] = e[1]
477+ */
478+ libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
479+ for (ii = 0; ii < 2; ii++) {
480+ w2[2][ii] = w2[3][ii];
481+ w2[1][ii] = w2[3][ii];
482+ }
483+ libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
484+ /*
485+ # IF(ANY(w2(1:2,1:4) < 0.0)):
486+ # WRITE(*,*) ie
487+ # WRITE(*,'(100e15.5)') x[0:4]
488+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
489+ # STOP "weighting 4=3=2"
490+ */
491+ }
492+ }
493+ else if (fabs(x[1] - x[0]) < thr) {
494+ /*
495+ e[3] = e[2], e[1] = e[0]
496+ */
497+ libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
498+ for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
499+ libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
500+ for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
501+ /*
502+ IF(ANY(w2(1:2,1:4) < 0.0)):
503+ WRITE(*,*) ie
504+ WRITE(*,'(100e15.5)') x[0:4]
505+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
506+ STOP "weighting 4=3 2=1"
507+ */
508+ }
509+ else {
510+ /*
511+ e[3] = e[2]
512+ */
513+ libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
514+ for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
515+ libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
516+ libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
517+ /*
518+ IF(ANY(w2(1:2,1:4) < 0.0)):
519+ WRITE(*,*) ie
520+ WRITE(*,'(100e15.5)') x[0:4]
521+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
522+ STOP "weighting 4=3"
523+ */
524+ }
525+ }
526+ else if (fabs(x[2] - x[1]) < thr) {
527+ if (fabs(x[2] - x[0]) < thr) {
528+ /*
529+ e[2] = e[1] = e[0]
530+ */
531+ libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
532+ libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
533+ for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
534+ for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[2][ii];
535+ /*
536+ IF(ANY(w2(1:2,1:4) < 0.0)):
537+ WRITE(*,*) ie
538+ WRITE(*,'(100e15.5)') x[0:4]
539+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
540+ STOP "weighting 3=2=1"
541+ */
542+ }
543+ else {
544+ /*
545+ e[2] = e[1]
546+ */
547+ libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
548+ libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
549+ for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
550+ libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
551+ /*
552+ IF(ANY(w2(1:2,1:4) < 0.0)):
553+ WRITE(*,*) ie
554+ WRITE(*,'(100e15.5)') x[0:4]
555+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
556+ STOP "weighting 3=2"
557+ */
558+ }
559+ }
560+ else if (fabs(x[1] - x[0]) < thr) {
561+ /*
562+ e[1] = e[0]
563+ */
564+ libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
565+ libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
566+ libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
567+ for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
568+ /*
569+ IF(ANY(w2(1:2,1:4) < 0.0)):
570+ WRITE(*,*) ie
571+ WRITE(*,'(100e15.5)') x[0:4]
572+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
573+ STOP "weighting 2=1"
574+ */
575+ }
576+ else {
577+ /*
578+ Different each other.
579+ */
580+ libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2], w2[3]);
581+ libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3], w2[2]);
582+ libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3], w2[1]);
583+ libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3], w2[0]);
584+ /*
585+ IF(ANY(w2(1:2,1:4) < 0.0)):
586+ WRITE(*,*) ie
587+ WRITE(*,'(100e15.5)') x[0:4]
588+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
589+ STOP "weighting"
590+ */
591+ }
592+ for (ii = 0; ii < 4; ii++) {
593+ w1[ie][indx[ii]][0] = w2[ii][0] / e0[ie][1];
594+ w1[ie][indx[ii]][1] = w2[ii][1] / (-e0[ie][1]);
595+ }
596+ }
597+}
598+//
599+// Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
600+// for 0<x<1, 0<y<1-x, 0<z<1-x-y
601+
602+
603+void libtetrabz_polcmplx_1234(
604+ double g1,
605+ double g2,
606+ double g3,
607+ double g4,
608+ double *w
609+ ) {
610+ /*
611+ 1, Different each other
612+ */
613+ double w2, w3, w4;
614+ /*
615+ Real
616+ */
617+ w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
618+ (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
619+ w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
620+ w2 = w2 / (g2 - g1);
621+ w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
622+ (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
623+ w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
624+ w3 = w3 / (g3 - g1);
625+ w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
626+ (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
627+ w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
628+ w4 = w4 / (g4 - g1);
629+ w2 = (w2 - w3) / (g2 - g3);
630+ w4 = (w4 - w3) / (g4 - g3);
631+ w[0] = (w4 - w2) / (2.0 * (g4 - g2));
632+ /*
633+ Imaginal
634+ */
635+ w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
636+ (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
637+ w2 = 4.0 * g2 - w2 / (g2 - g1);
638+ w2 = w2 / (g2 - g1);
639+ w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
640+ (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
641+ w3 = 4.0 * g3 - w3 / (g3 - g1);
642+ w3 = w3 / (g3 - g1);
643+ w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
644+ (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
645+ w4 = 4.0 * g4 - w4 / (g4 - g1);
646+ w4 = w4 / (g4 - g1);
647+ w2 = (w2 - w3) / (g2 - g3);
648+ w4 = (w4 - w3) / (g4 - g3);
649+ w[1] = (w4 - w2) / (2.0 * (g4 - g2));
650+}
651+
652+
653+void libtetrabz_polcmplx_1231(
654+ double g1,
655+ double g2,
656+ double g3,
657+ double *w
658+ ) {
659+ /*
660+ 2, g4 = g1
661+ */
662+ double w2, w3;
663+ /*
664+ Real
665+ */
666+ w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
667+ + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
668+ w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
669+ w2 = -g1 + w2 / (g2 - g1);
670+ w2 = w2 / (g2 - g1);
671+ w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
672+ + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
673+ w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
674+ w3 = -g1 + w3 / (g3 - g1);
675+ w3 = w3 / (g3 - g1);
676+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
677+ /*
678+ Imaginal
679+ */
680+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
681+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
682+ w2 = 4.0 * g2 - w2 / (g2 - g1);
683+ w2 = 1 + w2 / (g2 - g1);
684+ w2 = w2 / (g2 - g1);
685+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
686+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
687+ w3 = 4.0 * g3 - w3 / (g3 - g1);
688+ w3 = 1 + w3 / (g3 - g1);
689+ w3 = w3 / (g3 - g1);
690+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
691+}
692+
693+
694+void libtetrabz_polcmplx_1233(
695+ double g1,
696+ double g2,
697+ double g3,
698+ double *w
699+ ) {
700+ /*
701+ 3, g4 = g3
702+ */
703+ double w2, w3;
704+ /*
705+ Real
706+ */
707+ w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
708+ g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
709+ w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
710+ w2 = w2 / (g2 - g1);
711+ w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
712+ g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
713+ w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
714+ w3 = w3 / (g3 - g1);
715+ w2 = (w3 - w2) / (g3 - g2);
716+ w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
717+ + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3*g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
718+ w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
719+ w3 = 4.0 * g1 + w3 / (g3 - g1);
720+ w3 = w3 / (g3 - g1);
721+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
722+ /*
723+ Imaginal
724+ */
725+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
726+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
727+ w2 = 4.0 * g2 - w2 / (g2 - g1);
728+ w2 = w2 / (g2 - g1);
729+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
730+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
731+ w3 = 4.0 * g3 - w3 / (g3 - g1);
732+ w3 = w3 / (g3 - g1);
733+ w2 = (w3 - w2) / (g3 - g2);
734+ w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3*g3) * (atan(g3) - atan(g1))
735+ + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
736+ w3 = w3 / (g3 - g1) - 4.0 * g1;
737+ w3 = w3 / (g3 - g1) - 2.0;
738+ w3 = (2.0 * w3) / (g3 - g1);
739+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
740+}
741+
742+void libtetrabz_polcmplx_1221(
743+ double g1,
744+ double g2,
745+ double *w
746+ ) {
747+ /*
748+ 4, g4 = g1 and g3 = g2
749+ */
750+ /*
751+ Real
752+ */
753+ w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
754+ + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
755+ w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
756+ w[0] = 3.0 * g1 + w[0] / (g2 - g1);
757+ w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
758+ w[0] = w[0] / (2.0 * (g2 - g1));
759+ /*
760+ Imaginal
761+ */
762+ w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
763+ + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
764+ w[1] = -4.0 * g1 + w[1] / (g2 - g1);
765+ w[1] = -3.0 + w[1] / (g2 - g1);
766+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
767+}
768+
769+
770+void libtetrabz_polcmplx_1222(
771+ double g1,
772+ double g2,
773+ double *w
774+ ) {
775+ /*
776+ 5, g4 = g3 = g2
777+ */
778+ /*
779+ Real
780+ */
781+ w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
782+ + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
783+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
784+ w[0] = g1 - w[0] / (g2 - g1);
785+ w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
786+ w[0] = w[0] / (2.0 * (g2 - g1));
787+ /*
788+ Imaginal
789+ */
790+ w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
791+ + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
792+ w[1] = 4.0 * g1 + w[1] / (g2 - g1);
793+ w[1] = 1.0 + w[1] / (g2 - g1);
794+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
795+}
796+
797+
798+void libtetrabz_polcmplx_1211(
799+ double g1,
800+ double g2,
801+ double *w
802+) {
803+ /*
804+ 6, g4 = g3 = g1
805+ */
806+ /*
807+ Real
808+ */
809+ w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
810+ + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
811+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
812+ w[0] = -5.0 * g1 + w[0] / (g2 - g1);
813+ w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
814+ w[0] = w[0] / (6.0 * (g2 - g1));
815+ /*
816+ Imaginal
817+ */
818+ w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
819+ + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
820+ w[1] = 4.0 * g2 + w[1] / (g2 - g1);
821+ w[1] = 1.0 + w[1] / (g2 - g1);
822+ w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
823+}
--- /dev/null
+++ b/src/c/libtetrabz_polstat.c
@@ -0,0 +1,627 @@
1+//
2+// Copyright (c) 2014 Mitsuaki Kawamura
3+//
4+// Permission is hereby granted, free of charge, to any person obtaining a
5+// copy of this software and associated documentation files (the
6+// "Software"), to deal in the Software without restriction, including
7+// without limitation the rights to use, copy, modify, merge, publish,
8+// distribute, sublicense, and/or sell copies of the Software, and to
9+// permit persons to whom the Software is furnished to do so, subject to
10+// the following conditions:
11+//
12+// The above copyright notice and this permission notice shall be included
13+// in all copies or substantial portions of the Software.
14+//
15+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+//
23+#include "libtetrabz_common.h"
24+#include <math.h>
25+#include <stdio.h>
26+
27+/*
28+ Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
29+ for 0<x<1, 0<y<1-x, 0<z<1-x-y
30+*/
31+double libtetrabz_polstat_1234(
32+ double g1,
33+ double g2,
34+ double g3,
35+ double g4,
36+ double lng1,
37+ double lng2,
38+ double lng3,
39+ double lng4)
40+{
41+ /*
42+ 1, Different each other
43+ */
44+ double w2, w3, w4, w;
45+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 / (g2 - g1);
46+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 / (g3 - g1);
47+ w4 = ((lng4 - lng1) / (g4 - g1) * g4 - 1.0) * g4 / (g4 - g1);
48+ w2 = ((w2 - w3) * g2) / (g2 - g3);
49+ w4 = ((w4 - w3) * g4) / (g4 - g3);
50+ w = (w4 - w2) / (g4 - g2);
51+ return w;
52+}
53+
54+
55+double libtetrabz_polstat_1231(
56+ double g1,
57+ double g2,
58+ double g3,
59+ double lng1,
60+ double lng2,
61+ double lng3
62+) {
63+ /*
64+ 2, g4 = g1
65+ */
66+ double w2, w3, w;
67+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 * g2 / (g2 - g1) - g1 / 2.0;
68+ w2 = w2 / (g2 - g1);
69+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 * g3 / (g3 - g1) - g1 / 2.0;
70+ w3 = w3 / (g3 - g1);
71+ w = (w3 - w2) / (g3 - g2);
72+ return w;
73+}
74+
75+double libtetrabz_polstat_1233(
76+ double g1,
77+ double g2,
78+ double g3,
79+ double lng1,
80+ double lng2,
81+ double lng3
82+) {
83+ /*
84+ # 3, g4 = g3
85+ */
86+ double w2, w3, w;
87+ w2 = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
88+ w2 = (g2 * w2) / (g2 - g1);
89+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
90+ w3 = (g3 * w3) / (g3 - g1);
91+ w2 = (w3 - w2) / (g3 - g2);
92+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
93+ w3 = 1.0 - (2.0 * w3 * g1) / (g3 - g1);
94+ w3 = w3 / (g3 - g1);
95+ w = (g3 * w3 - g2 * w2) / (g3 - g2);
96+ return w;
97+}
98+
99+double libtetrabz_polstat_1221(
100+ double g1,
101+ double g2,
102+ double lng1,
103+ double lng2
104+) {
105+ /*
106+ 4, g4 = g1 and g3 = g2
107+ */
108+ double w;
109+ w = 1.0 - (lng2 - lng1) / (g2 - g1) * g1;
110+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
111+ w = -1.0 + (3.0 * g2 * w) / (g2 - g1);
112+ w = w / (2.0 * (g2 - g1));
113+ return w;
114+}
115+
116+double libtetrabz_polstat_1222(
117+ double g1,
118+ double g2,
119+ double lng1,
120+ double lng2
121+) {
122+ /*
123+ 5, g4 = g3 = g2
124+ */
125+ double w;
126+ w = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
127+ w = (2.0 * g1 * w) / (g2 - g1) - 1.0;
128+ w = (3.0 * g1 * w) / (g2 - g1) + 1.0;
129+ w = w / (2.0 * (g2 - g1));
130+ return w;
131+}
132+
133+double libtetrabz_polstat_1211(
134+ double g1,
135+ double g2,
136+ double lng1,
137+ double lng2
138+) {
139+ /*
140+ 6, g4 = g3 = g1
141+ */
142+ double w;
143+ w = -1.0 + (lng2 - lng1) / (g2 - g1) * g2;
144+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
145+ w = -1.0 + (3.0 * g2 * w) / (2.0 * (g2 - g1));
146+ w = w / (3.0 * (g2 - g1));
147+ return w;
148+}
149+
150+void libtetrabz_polstat3(
151+ double *de,
152+ double *w1
153+)
154+{
155+ /*
156+ Tetrahedron method for delta(om - ep + e)
157+ */
158+ int ii, indx[4];
159+ double thr, thr2, e[4], ln[4];
160+
161+ for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
162+ eig_sort(4, e, indx);
163+
164+ thr = e[3] * 1.0e-3;
165+ thr2 = 1.0e-8;
166+
167+ for(ii=0;ii<4;ii++){
168+ if (e[ii] < thr2) {
169+ if (ii == 3) {
170+ printf(" Nesting # ");
171+ }
172+ ln[ii] = 0.0;
173+ e[ii] = 0.0;
174+ }
175+ else{
176+ ln[ii] = log(e[ii]);
177+ }
178+ }
179+
180+ if (fabs(e[3] - e[2]) < thr) {
181+ if (fabs(e[3] - e[1]) < thr) {
182+ if (fabs(e[3] - e[0]) < thr) {
183+ /*
184+ e[3] = e[2] = e[1] = e[0]
185+ */
186+ w1[indx[3]] = 0.25 / e[3];
187+ w1[indx[2]] = w1[indx[3]];
188+ w1[indx[1]] = w1[indx[3]];
189+ w1[indx[0]] = w1[indx[3]];
190+ }
191+ else {
192+ /*
193+ e[3] = e[2] = e[1]
194+ */
195+ w1[indx[3]] = libtetrabz_polstat_1211(e[3], e[0], ln[3], ln[0]);
196+ w1[indx[2]] = w1[indx[3]];
197+ w1[indx[1]] = w1[indx[3]];
198+ w1[indx[0]] = libtetrabz_polstat_1222(e[0], e[3], ln[0], ln[3]);
199+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
200+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
201+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
202+ printf("weighting 4=3=2");
203+ }
204+ }
205+ }
206+ else if (fabs(e[1] - e[0]) < thr) {
207+ /*
208+ e[3] = e[2], e[1] = e[0]
209+ */
210+ w1[indx[3]] = libtetrabz_polstat_1221(e[3], e[1], ln[3], ln[1]);
211+ w1[indx[2]] = w1[indx[3]];
212+ w1[indx[1]] = libtetrabz_polstat_1221(e[1], e[3], ln[1], ln[3]);
213+ w1[indx[0]] = w1[indx[1]];
214+
215+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
216+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
217+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
218+ printf("weighting 4=3 2=1");
219+ }
220+ }
221+ else {
222+ /*
223+ e[3] = e[2]
224+ */
225+ w1[indx[3]] = libtetrabz_polstat_1231(e[3], e[0], e[1], ln[3], ln[0], ln[1]);
226+ w1[indx[2]] = w1[indx[3]];
227+ w1[indx[1]] = libtetrabz_polstat_1233(e[1], e[0], e[3], ln[1], ln[0], ln[3]);
228+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[1], e[3], ln[0], ln[1], ln[3]);
229+
230+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
231+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
232+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
233+ printf("weighting 4=3");
234+ }
235+ }
236+ }
237+ else if(fabs(e[2] - e[1]) < thr){
238+ if (fabs(e[2] - e[0]) < thr) {
239+ /*
240+ e[2] = e[1] = e[0]
241+ */
242+ w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
243+ w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
244+ w1[indx[1]] = w1[indx[2]];
245+ w1[indx[0]] = w1[indx[2]];
246+
247+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
248+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
249+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
250+ printf("weighting 3=2=1");
251+ }
252+ }
253+ else {
254+ /*
255+ e[2] = e[1]
256+ */
257+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
258+ w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
259+ w1[indx[1]] = w1[indx[2]];
260+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
261+
262+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
263+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
264+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
265+ printf("weighting 3=2");
266+ }
267+ }
268+ }
269+ else if(fabs(e[1] - e[0]) < thr) {
270+ /*
271+ e[1] = e[0]
272+ */
273+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
274+ w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
275+ w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
276+ w1[indx[0]] = w1[indx[1]];
277+
278+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
279+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
280+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
281+ printf("weighting 2=1");
282+ }
283+ }
284+ else {
285+ /*
286+ Different each other.
287+ */
288+ w1[indx[3]] = libtetrabz_polstat_1234(e[3], e[0], e[1], e[2], ln[3], ln[0], ln[1], ln[2]);
289+ w1[indx[2]] = libtetrabz_polstat_1234(e[2], e[0], e[1], e[3], ln[2], ln[0], ln[1], ln[3]);
290+ w1[indx[1]] = libtetrabz_polstat_1234(e[1], e[0], e[2], e[3], ln[1], ln[0], ln[2], ln[3]);
291+ w1[indx[0]] = libtetrabz_polstat_1234(e[0], e[1], e[2], e[3], ln[0], ln[1], ln[2], ln[3]);
292+
293+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
294+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
295+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
296+ printf("weighting");
297+ }
298+ }
299+}
300+
301+void libtetrabz_polstat2(
302+ int nb,
303+ double *ei1,
304+ double **ej1,
305+ double **w1
306+) {
307+ /*
308+ Tetrahedron method for theta( - E2)
309+ :param ei1:
310+ :param ej1:
311+ :return:
312+ */
313+ int ii, jj, ib, indx[4];
314+ double v, thr, e[4], de[4], w2[4], tsmall[4][4];
315+
316+ thr = 1.0e-8;
317+
318+ for (ib = 0; ib < nb; ib++) {
319+
320+ for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
321+ eig_sort(4, e, indx);
322+
323+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
324+
325+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
326+
327+ if (v > thr) {
328+ for (ii = 0; ii < 4; ii++) {
329+ de[ii] = 0.0;
330+ for (jj = 0; jj < 4; jj++)
331+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
332+ }
333+ libtetrabz_polstat3(de, w2);
334+ for (ii = 0; ii < 4; ii++)
335+ for (jj = 0; jj < 4; jj++)
336+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
337+ }
338+ }
339+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
340+
341+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
342+
343+ if (v > thr) {
344+ for (ii = 0; ii < 4; ii++) {
345+ de[ii] = 0.0;
346+ for (jj = 0; jj < 4; jj++)
347+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
348+ }
349+ libtetrabz_polstat3(de, w2);
350+ for (ii = 0; ii < 4; ii++)
351+ for (jj = 0; jj < 4; jj++)
352+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
353+ }
354+
355+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
356+
357+ if (v > thr) {
358+ for (ii = 0; ii < 4; ii++) {
359+ de[ii] = 0.0;
360+ for (jj = 0; jj < 4; jj++)
361+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
362+ }
363+ libtetrabz_polstat3(de, w2);
364+ for (ii = 0; ii < 4; ii++)
365+ for (jj = 0; jj < 4; jj++)
366+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
367+ }
368+
369+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
370+
371+ if (v > thr) {
372+ for (ii = 0; ii < 4; ii++) {
373+ de[ii] = 0.0;
374+ for (jj = 0; jj < 4; jj++)
375+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
376+ }
377+ libtetrabz_polstat3(de, w2);
378+ for (ii = 0; ii < 4; ii++)
379+ for (jj = 0; jj < 4; jj++)
380+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
381+ }
382+ }
383+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
384+
385+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
386+
387+ if (v > thr) {
388+ for (ii = 0; ii < 4; ii++) {
389+ de[ii] = 0.0;
390+ for (jj = 0; jj < 4; jj++)
391+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
392+ }
393+ libtetrabz_polstat3(de, w2);
394+ for (ii = 0; ii < 4; ii++)
395+ for (jj = 0; jj < 4; jj++)
396+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
397+ }
398+
399+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
400+
401+ if (v > thr) {
402+ for (ii = 0; ii < 4; ii++) {
403+ de[ii] = 0.0;
404+ for (jj = 0; jj < 4; jj++)
405+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
406+ }
407+ libtetrabz_polstat3(de, w2);
408+ for (ii = 0; ii < 4; ii++)
409+ for (jj = 0; jj < 4; jj++)
410+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
411+ }
412+
413+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
414+
415+ if (v > thr) {
416+ for (ii = 0; ii < 4; ii++) {
417+ de[ii] = 0.0;
418+ for (jj = 0; jj < 4; jj++)
419+ de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
420+ }
421+ libtetrabz_polstat3(de, w2);
422+ for (ii = 0; ii < 4; ii++)
423+ for (jj = 0; jj < 4; jj++)
424+ w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
425+ }
426+ }
427+ else if (e[3] <= 0.0) {
428+ for (ii = 0; ii < 4; ii++)
429+ de[ii] = ej1[indx[ii]][ib] - ei1[indx[ii]];
430+ libtetrabz_polstat3(de, w2);
431+ for (ii = 0; ii < 4; ii++) w1[ib][ii] += w2[ii];
432+ }
433+ }
434+}
435+
436+void polstat(
437+ int *ng,
438+ int nk,
439+ int nb,
440+ double **bvec,
441+ double **eig1,
442+ double **eig2,
443+ double ***wght
444+)
445+{
446+ /*
447+ Compute Static polarization function
448+ */
449+ int it, ik, ib, ii, jj, jb, **ikv, indx[4];
450+ double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
451+
452+ thr = 1.0e-10;
453+
454+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
455+
456+ for (ik = 0; ik < nk; ik++)
457+ for (ib = 0; ib < nb; ib++)
458+ for (jb = 0; jb < nb; jb++)
459+ wght[ik][ib][jb] = 0.0;
460+
461+ for(it = 0; it < 6*nk; it++) {
462+
463+ for (ii = 0; ii < 4; ii++) {
464+ for (ib = 0; ib < nb; ib++) {
465+ ei1[ii][ib] = 0.0;
466+ ej1[ii][ib] = 0.0;
467+ for (jj = 0; jj < 20; jj++) {
468+ ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
469+ ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
470+ }
471+ }
472+ }
473+
474+ for (ib = 0; ib < nb; ib++)
475+ for (jb = 0; jb < nb; jb++)
476+ for (ii = 0; ii < 4; ii++)
477+ w1[ib][jb][ii] = 0.0;
478+
479+ for (ib = 0; ib < nb; ib++) {
480+
481+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
482+ eig_sort(4, e, indx);
483+
484+ if (e[0] <= 0.0 && 0.0 < e[1]) {
485+
486+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
487+
488+ if (v > thr) {
489+ for (ii = 0; ii < 4; ii++)
490+ for (jj = 0; jj < 4; jj++) {
491+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
492+ for (jb = 0; jb < nb; jb++)
493+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
494+ }
495+ libtetrabz_polstat2(nb, ei2, ej2, w2);
496+ for (ii = 0; ii < 4; ii++)
497+ for (jj = 0; jj < 4; jj++)
498+ for (jb = 0; jb < nb; jb++)
499+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
500+ }
501+ }
502+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
503+
504+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
505+
506+ if (v > thr) {
507+ for (ii = 0; ii < 4; ii++)
508+ for (jj = 0; jj < 4; jj++) {
509+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
510+ for (jb = 0; jb < nb; jb++)
511+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
512+ }
513+ libtetrabz_polstat2(nb, ei2, ej2, w2);
514+ for (ii = 0; ii < 4; ii++)
515+ for (jj = 0; jj < 4; jj++)
516+ for (jb = 0; jb < nb; jb++)
517+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
518+ }
519+
520+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
521+
522+ if (v > thr) {
523+ for (ii = 0; ii < 4; ii++)
524+ for (jj = 0; jj < 4; jj++) {
525+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
526+ for (jb = 0; jb < nb; jb++)
527+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
528+ }
529+ libtetrabz_polstat2(nb, ei2, ej2, w2);
530+ for (ii = 0; ii < 4; ii++)
531+ for (jj = 0; jj < 4; jj++)
532+ for (jb = 0; jb < nb; jb++)
533+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
534+ }
535+
536+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
537+
538+ if (v > thr) {
539+ for (ii = 0; ii < 4; ii++)
540+ for (jj = 0; jj < 4; jj++) {
541+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
542+ for (jb = 0; jb < nb; jb++)
543+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
544+ }
545+ libtetrabz_polstat2(nb, ei2, ej2, w2);
546+ for (ii = 0; ii < 4; ii++)
547+ for (jj = 0; jj < 4; jj++)
548+ for (jb = 0; jb < nb; jb++)
549+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
550+ }
551+ }
552+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
553+
554+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
555+
556+ if (v > thr) {
557+ for (ii = 0; ii < 4; ii++)
558+ for (jj = 0; jj < 4; jj++) {
559+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
560+ for (jb = 0; jb < nb; jb++)
561+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
562+ }
563+ libtetrabz_polstat2(nb, ei2, ej2, w2);
564+ for (ii = 0; ii < 4; ii++)
565+ for (jj = 0; jj < 4; jj++)
566+ for (jb = 0; jb < nb; jb++)
567+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
568+ }
569+
570+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
571+
572+ if (v > thr) {
573+ for (ii = 0; ii < 4; ii++)
574+ for (jj = 0; jj < 4; jj++) {
575+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
576+ for (jb = 0; jb < nb; jb++)
577+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
578+ }
579+ libtetrabz_polstat2(nb, ei2, ej2, w2);
580+ for (ii = 0; ii < 4; ii++)
581+ for (jj = 0; jj < 4; jj++)
582+ for (jb = 0; jb < nb; jb++)
583+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
584+ }
585+
586+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
587+
588+ if (v > thr) {
589+ for (ii = 0; ii < 4; ii++)
590+ for (jj = 0; jj < 4; jj++) {
591+ ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
592+ for (jb = 0; jb < nb; jb++)
593+ ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
594+ }
595+ libtetrabz_polstat2(nb, ei2, ej2, w2);
596+ for (ii = 0; ii < 4; ii++)
597+ for (jj = 0; jj < 4; jj++)
598+ for (jb = 0; jb < nb; jb++)
599+ w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
600+ }
601+ }
602+ else if (e[3] <= 0.0) {
603+ for (ii = 0; ii < 4; ii++) {
604+ ei2[ii] = ej1[ii][ib];
605+ for (jb = 0; jb < nb; jb++)
606+ ej2[ii][jb] = ej1[ii][jb];
607+ }
608+ libtetrabz_polstat2(nb, ei2, ej2, w2);
609+ for (ii = 0; ii < 4; ii++)
610+ for (jb = 0; jb < nb; jb++)
611+ w1[ib][jb][ii] += v * w2[jb][ii];
612+ }
613+ else {
614+ continue;
615+ }
616+ }
617+ for (ii = 0; ii < 20; ii++)
618+ for (jj = 0; jj < 4; jj++)
619+ for (ib = 0; ib < nb; ib++)
620+ for (jb = 0; jb < nb; jb++)
621+ wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
622+ }
623+ for (ik = 0; ik < nk; ik++)
624+ for (ib = 0; ib < nb; ib++)
625+ for (jb = 0; jb < nb; jb++)
626+ wght[ik][ib][jb] /= (6.0 * (double) nk);
627+}
Show on old repository browser