OSDN -- オープンソース・ソフトウェアの開発とダウンロード
日本語の翻訳状況 translation status for ja
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

python: コミット

libtetrabz python package


コミットメタ情報

リビジョンe929f2556ff90ade558a10373b57a9642ead72bc (tree)
日時2022-03-24 20:54:49
作者Mitsuaki Kawamura <kawamitsuaki@gmai...>
コミッターMitsuaki Kawamura

ログメッセージ

These files are integrated into libtetrabz.c

変更サマリ

  • delete: src/c/__init__.c
  • delete: src/c/libtetrabz_common.c
  • delete: src/c/libtetrabz_common.h
  • delete: src/c/libtetrabz_dbldelta.c
  • delete: src/c/libtetrabz_dblstep.c
  • delete: src/c/libtetrabz_dos.c
  • delete: src/c/libtetrabz_fermigr.c
  • delete: src/c/libtetrabz_occ.c
  • delete: src/c/libtetrabz_polcmplx.c
  • delete: src/c/libtetrabz_polstat.c
  • delete: src/c/test.c
  • delete: src/c/test.py

差分

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