• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revision9ea7b0727c5abc70b55c6026f6faee484d624519 (tree)
Zeit2021-11-05 00:49:52
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Add test which is the same as the fortran ver.

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/src/dos.py
@@ -0,0 +1,32 @@
1+import time
2+
3+import numpy
4+import libtetrabz
5+import matplotlib.pyplot as plt
6+import matplotlib
7+
8+
9+bvec = numpy.array([[10.0, 0.0, 0.0],
10+ [0.0, 10.0, 0.0],
11+ [0.0, 0.0, 10.0]])
12+
13+ng0 = 10
14+nb = 1
15+ng = numpy.array([ng0, ng0, ng0])
16+eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
17+for i0 in range(ng[0]):
18+ for i1 in range(ng[1]):
19+ for i2 in range(ng[2]):
20+ kvec = 2.0*numpy.pi*numpy.array([i0, i1, i2]) / ng[0:3]
21+ eig[i0, i1, i2, 0] = - numpy.cos(kvec).sum(0)
22+
23+e0 = numpy.linspace(-3, 3, 100)
24+
25+t0 = time.time()
26+wght = libtetrabz.dos(bvec, eig, e0)
27+t1 = time.time()
28+print("Time", t1-t0, "sec", wght.sum())
29+dos = wght.sum(3).sum(2).sum(1).sum(0)
30+matplotlib.use('TkAgg')
31+plt.plot(e0, dos)
32+plt.show()
--- a/src/libtetrabz/libtetrabz_dbldelta.py
+++ b/src/libtetrabz/libtetrabz_dbldelta.py
@@ -51,9 +51,10 @@ def dbldelta(bvec, eig1, eig2):
5151 ej1 = wlsm.dot(eig2t)
5252 #
5353 w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
54+ w2t = numpy.zeros([nb, 4], dtype=numpy.float_)
5455 for ib in range(nb):
5556 #
56- e = ei1[0:4, ib]
57+ e = ei1[0:4, ib].copy()
5758 indx = e.argsort(0)
5859 e.sort(0)
5960 #
@@ -62,10 +63,10 @@ def dbldelta(bvec, eig1, eig2):
6263 v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e)
6364 #
6465 if v > thr:
65- #
6666 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
6767 w2 = libtetrabz_dbldelta2(ej2)
68- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
68+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
69+ w1[ib, 0:nb, 0:4] += v * w2t
6970 #
7071 elif e[1] < 0.0 <= e[2]:
7172 #
@@ -74,14 +75,16 @@ def dbldelta(bvec, eig1, eig2):
7475 if v > thr:
7576 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7677 w2 = libtetrabz_dbldelta2(ej2)
77- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
78+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
79+ w1[ib, 0:nb, 0:4] += v * w2t
7880 #
7981 v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e)
8082 #
8183 if v > thr:
8284 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8385 w2 = libtetrabz_dbldelta2(ej2)
84- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
86+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
87+ w1[ib, 0:nb, 0:4] += v * w2t
8588 #
8689 elif e[2] < 0.0 < e[3]:
8790 #
@@ -90,7 +93,8 @@ def dbldelta(bvec, eig1, eig2):
9093 if v > thr:
9194 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9295 w2 = libtetrabz_dbldelta2(ej2)
93- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
96+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
97+ w1[ib, 0:nb, 0:4] += v * w2t
9498 #
9599 else:
96100 continue
@@ -110,35 +114,40 @@ def libtetrabz_dbldelta2(ej):
110114 """
111115 nb = ej.shape[1]
112116 w = numpy.zeros([nb, 3], dtype=numpy.float_)
113- a = numpy.empty([3, 3], dtype=numpy.float_)
114117 for ib in range(nb):
115118 #
116119 if abs(ej[0:3, ib]).max() < 1.0e-10:
117120 raise ValueError("Nesting ##")
118121 #
119- e = ej[0:3, ib]
122+ e = ej[0:3, ib].copy()
120123 indx = e.argsort(0)
121124 e.sort(0)
122125 #
123- for ii in range(3):
124- a[0:3, ii] = (0.0 - e[ii]) / (e[0:3] - e[ii])
126+ #for ii in range(3):
127+ # a[0:3, ii] = (0.0 - e[ii]) / (e[0:3] - e[ii])
125128 #
126129 if e[0] < 0.0 <= e[1] or e[0] <= 0.0 < e[1]:
127130 #
131+ a10 = (0.0 - e[0]) / (e[1] - e[0])
132+ a20 = (0.0 - e[0]) / (e[2] - e[0])
133+ #
128134 # v = a[1, 0] * a[2, 0] / (0.0 - e[0])
129- v = a[1, 0] / (e[2] - e[0])
135+ v = a10 / (e[2] - e[0])
130136 #
131- w[ib, indx[0]] = v * (a[0, 1] + a[0, 2])
132- w[ib, indx[1]] = v * a[1, 0]
133- w[ib, indx[2]] = v * a[2, 0]
137+ w[ib, indx[0]] = v * (2.0 - a10 - a20)
138+ w[ib, indx[1]] = v * a10
139+ w[ib, indx[2]] = v * a20
134140 #
135141 elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
136142 #
143+ a02 = (0.0 - e[2]) / (e[0] - e[2])
144+ a12 = (0.0 - e[2]) / (e[1] - e[2])
145+ #
137146 # v = a[1, 2] * a[1, 2] / (e[2] - 0.0)
138- v = a[1, 2] / (e[2] - e[0])
147+ v = a12 / (e[2] - e[0])
139148 #
140- w[ib, indx[0]] = v * a[0, 2]
141- w[ib, indx[1]] = v * a[1, 2]
142- w[ib, indx[2]] = v * (a[2, 0] + a[2, 1])
149+ w[ib, indx[0]] = v * a02
150+ w[ib, indx[1]] = v * a12
151+ w[ib, indx[2]] = v * (2.0 - a02 - a12)
143152 #
144153 return w
--- a/src/libtetrabz/libtetrabz_dblstep.py
+++ b/src/libtetrabz/libtetrabz_dblstep.py
@@ -52,13 +52,14 @@ def dblstep(bvec, eig1, eig2):
5252 ej1 = wlsm.dot(eig2t)
5353 #
5454 w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
55+ w2t = numpy.zeros([nb, 4], dtype=numpy.float_)
5556 for ib in range(nb):
5657 #
57- e = ei1[0:4, ib]
58+ e = ei1[0:4, ib].copy()
5859 indx = e.argsort(0)
5960 e.sort(0)
6061 #
61- if e(1) <= 0.0 < e(2):
62+ if e[0] <= 0.0 < e[1]:
6263 #
6364 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
6465 #
@@ -66,9 +67,10 @@ def dblstep(bvec, eig1, eig2):
6667 ei2 = tsmall.dot(ei1[indx[0:4], ib])
6768 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
6869 w2 = libtetrabz_dblstep2(ei2, ej2)
69- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
70+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
71+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
7072 #
71- elif e(2) <= 0.0 < e(3):
73+ elif e[1] <= 0.0 < e[2]:
7274 #
7375 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
7476 #
@@ -76,7 +78,8 @@ def dblstep(bvec, eig1, eig2):
7678 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7779 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7880 w2 = libtetrabz_dblstep2(ei2, ej2)
79- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
81+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
82+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
8083 #
8184 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8285 #
@@ -84,7 +87,8 @@ def dblstep(bvec, eig1, eig2):
8487 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8588 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8689 w2 = libtetrabz_dblstep2(ei2, ej2)
87- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
90+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
91+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
8892 #
8993 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9094 #
@@ -92,9 +96,10 @@ def dblstep(bvec, eig1, eig2):
9296 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9397 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9498 w2 = libtetrabz_dblstep2(ei2, ej2)
95- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
99+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
100+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
96101 #
97- elif e(3) <= 0.0 < e(4):
102+ elif e[2] <= 0.0 < e[3]:
98103 #
99104 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
100105 #
@@ -102,7 +107,8 @@ def dblstep(bvec, eig1, eig2):
102107 ei2 = tsmall.dot(ei1[indx[0:4], ib])
103108 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
104109 w2 = libtetrabz_dblstep2(ei2, ej2)
105- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
110+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
111+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
106112 #
107113 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
108114 #
@@ -110,7 +116,8 @@ def dblstep(bvec, eig1, eig2):
110116 ei2 = tsmall.dot(ei1[indx[0:4], ib])
111117 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
112118 w2 = libtetrabz_dblstep2(ei2, ej2)
113- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
119+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
120+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
114121 #
115122 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
116123 #
@@ -118,9 +125,10 @@ def dblstep(bvec, eig1, eig2):
118125 ei2 = tsmall.dot(ei1[indx[0:4], ib])
119126 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
120127 w2 = libtetrabz_dblstep2(ei2, ej2)
121- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
128+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
129+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
122130 #
123- elif e(4) <= 0.0:
131+ elif e[3] <= 0.0:
124132 ei2 = ei1[0:4, ib]
125133 ej2 = ej1[0:4, 0:nb]
126134 w2 = libtetrabz_dblstep2(ei2, ej2)
@@ -152,19 +160,19 @@ def libtetrabz_dblstep2(ei1, ej1):
152160 indx = e.argsort(0)
153161 e.sort(0)
154162 #
155- if abs(e(1)) < thr and abs(e(4)) < thr:
163+ if abs(e[0]) < thr and abs(e[3]) < thr:
156164 #
157165 # Theta(0) = 0.5
158166 #
159167 v = 0.5
160168 w1[ib, 0:4] += v * 0.25
161169 #
162- elif e(1) <= 0.0 < e(2) or e(1) < 0.0 <= e(2):
170+ elif e[0] <= 0.0 < e[1] or e[0] < 0.0 <= e[1]:
163171 #
164172 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
165173 w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
166174 #
167- elif e(2) <= 0.0 < e(3) or e(2) < 0.0 <= e(3):
175+ elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
168176 #
169177 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
170178 w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
@@ -175,7 +183,7 @@ def libtetrabz_dblstep2(ei1, ej1):
175183 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
176184 w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
177185 #
178- elif e(3) <= 0.0 < e(4) or e(3) < 0.0 <= e(4):
186+ elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
179187 #
180188 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
181189 w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
@@ -186,7 +194,7 @@ def libtetrabz_dblstep2(ei1, ej1):
186194 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
187195 w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
188196 #
189- elif e(4) <= 0.0:
197+ elif e[3] <= 0.0:
190198 #
191199 w1[ib, 0:4] += 0.25
192200 #
--- a/src/libtetrabz/libtetrabz_dos.py
+++ b/src/libtetrabz/libtetrabz_dos.py
@@ -32,15 +32,14 @@ def dos(bvec, eig, e0):
3232 :param ndarray(ne, float) e0: Energy grid
3333 :return ndarray([ng0, ng1, ng2, nb, ne], float) wght: Weight
3434 """
35-
3635 ng = numpy.array(eig.shape[0:3])
3736 nk = ng.prod(0)
3837 nb = eig.shape[3]
3938 ne = e0.shape[0]
4039 wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
41-
40+ #
4241 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, ne], dtype=numpy.float_)
43-
42+ #
4443 for it in range(6*nk):
4544 #
4645 eigt = numpy.empty([20, nb], dtype=numpy.float_)
@@ -52,7 +51,7 @@ def dos(bvec, eig, e0):
5251 #
5352 for ib in range(nb):
5453 #
55- e = ei1[0:4, ib]
54+ e = ei1[0:4, ib].copy()
5655 indx = e.argsort(0)
5756 e = e[indx[0:4]]
5857 #
@@ -116,7 +115,7 @@ def intdos(bvec, eig, e0):
116115 #
117116 for ib in range(nb):
118117 #
119- e = ei1[0:4, ib]
118+ e = ei1[0:4, ib].copy()
120119 indx = e.argsort(0)
121120 e.sort(0)
122121 #
@@ -138,7 +137,7 @@ def intdos(bvec, eig, e0):
138137 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e - e0[ie])
139138 w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
140139 #
141- elif e[2] <= e0[ie] < e[3] or e[2] < e0[ie] .AND. e0[ie] <= e[3]:
140+ elif e[2] <= e0[ie] < e[3] or e[2] < e0[ie] <= e[3]:
142141 #
143142 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e - e0[ie])
144143 w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
@@ -149,7 +148,7 @@ def intdos(bvec, eig, e0):
149148 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e - e0[ie])
150149 w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
151150 #
152- elif e[3] <= e0(ie):
151+ elif e[3] <= e0[ie]:
153152 w1[ib, ie, 0:4] = 0.25
154153 else:
155154 continue
--- a/src/libtetrabz/libtetrabz_fermigr.py
+++ b/src/libtetrabz/libtetrabz_fermigr.py
@@ -54,9 +54,10 @@ def fermigr(bvec, eig1, eig2, e0):
5454 ej1 = wlsm.dot(eig2t)
5555 #
5656 w1 = numpy.zeros([nb, nb, ne, 4], dtype=numpy.float_)
57+ w2t = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
5758 for ib in range(nb):
5859 #
59- e = ei1[0:4, ib]
60+ e = ei1[0:4, ib].copy()
6061 indx = e.argsort(0)
6162 e.sort(0)
6263 #
@@ -69,7 +70,8 @@ def fermigr(bvec, eig1, eig2, e0):
6970 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7071 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7172 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
72- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
73+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
74+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
7375 #
7476 elif e[1] <= 0.0 < e[2]:
7577 #
@@ -80,7 +82,8 @@ def fermigr(bvec, eig1, eig2, e0):
8082 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8183 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8284 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
83- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
85+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
86+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
8487 #
8588 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8689 #
@@ -89,7 +92,8 @@ def fermigr(bvec, eig1, eig2, e0):
8992 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9093 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9194 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
92- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
95+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
96+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
9397 #
9498 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9599 #
@@ -98,7 +102,8 @@ def fermigr(bvec, eig1, eig2, e0):
98102 ei2 = tsmall.dot(ei1[indx[0:4], ib])
99103 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
100104 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
101- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
105+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
106+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
102107 #
103108 elif e[2] <= 0.0 < e[3]:
104109 #
@@ -109,7 +114,8 @@ def fermigr(bvec, eig1, eig2, e0):
109114 ei2 = tsmall.dot(ei1[indx[0:4], ib])
110115 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
111116 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
112- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
117+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
118+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
113119 #
114120 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
115121 #
@@ -118,7 +124,8 @@ def fermigr(bvec, eig1, eig2, e0):
118124 ei2 = tsmall.dot(ei1[indx[0:4], ib])
119125 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
120126 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
121- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
127+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
128+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
122129 #
123130 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
124131 #
@@ -127,7 +134,8 @@ def fermigr(bvec, eig1, eig2, e0):
127134 ei2 = tsmall.dot(ei1[indx[0:4], ib])
128135 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
129136 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
130- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
137+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
138+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
131139 #
132140 elif e[3] <= 0.0:
133141 #
@@ -158,6 +166,7 @@ def libtetrabz_fermigr2(e0, ei1, ej1):
158166 ne = e0.shape[0]
159167 nb = ej1.shape[1]
160168 w1 = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
169+ w2t = numpy.zeros([ne, 4], dtype=numpy.float_)
161170 thr = 1.0e-8
162171 #
163172 for ib in range(nb):
@@ -174,7 +183,8 @@ def libtetrabz_fermigr2(e0, ei1, ej1):
174183 #
175184 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
176185 w2 = libtetrabz_fermigr3(e0, de)
177- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
186+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
187+ w1[ib, 0:ne, 0:4] += v * w2t
178188 #
179189 #
180190 elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
@@ -182,61 +192,54 @@ def libtetrabz_fermigr2(e0, ei1, ej1):
182192 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
183193 #
184194 if v > thr:
185- #
186195 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
187196 w2 = libtetrabz_fermigr3(e0, de)
188- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
189- #
197+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
198+ w1[ib, 0:ne, 0:4] += v * w2t
190199 #
191200 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
192201 #
193202 if v > thr:
194- #
195203 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
196204 w2 = libtetrabz_fermigr3(e0, de)
197- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
198- #
205+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
206+ w1[ib, 0:ne, 0:4] += v * w2t
199207 #
200208 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
201209 #
202210 if v > thr:
203- #
204211 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
205212 w2 = libtetrabz_fermigr3(e0, de)
206- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
207- #
213+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
214+ w1[ib, 0:ne, 0:4] += v * w2t
208215 #
209216 elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
210217 #
211218 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
212219 #
213220 if v > thr:
214- #
215221 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
216222 w2 = libtetrabz_fermigr3(e0, de)
217- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
218- #
223+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
224+ w1[ib, 0:ne, 0:4] += v * w2t
219225 #
220226 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
221227 #
222228 if v > thr:
223- #
224229 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
225230 w2 = libtetrabz_fermigr3(e0, de)
226- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
227- #
231+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
232+ w1[ib, 0:ne, 0:4] += v * w2t
228233 #
229234 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
230235 #
231236 if v > thr:
232- #
233237 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
234238 w2 = libtetrabz_fermigr3(e0, de)
235- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
236- #
239+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
240+ w1[ib, 0:ne, 0:4] += v * w2t
237241 #
238242 elif e[3] <= 0.0:
239- #
240243 de = ej1[0:4, ib] - ei1[0:4]
241244 w2 = libtetrabz_fermigr3(e0, de)
242245 w1[ib, 0:ne, 0:4] += w2
@@ -247,19 +250,19 @@ def libtetrabz_fermigr2(e0, ei1, ej1):
247250 def libtetrabz_fermigr3(e0, de):
248251 #
249252 ne = e0.shape[0]
250- e = de[0:4]
253+ e = de[0:4].copy()
251254 indx = e.argsort(0)
252255 e.sort(0)
253256 w1 = numpy.zeros([ne, 4], dtype=numpy.float_)
254257 #
255258 for ie in range(ne):
256259 #
257- if e[0] < e0[ie] .AND. e0[ie] <= e[1]:
260+ if e[0] < e0[ie] <= e[1]:
258261 #
259262 v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e[0:4] - e0[ie])
260263 w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
261264 #
262- elif e[1] < e0[ie] .AND. e0[ie] <= e[2]:
265+ elif e[1] < e0[ie] <= e[2]:
263266 #
264267 v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e[0:4] - e0[ie])
265268 w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
@@ -267,7 +270,7 @@ def libtetrabz_fermigr3(e0, de):
267270 v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e[0:4] - e0[ie])
268271 w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
269272 #
270- elif e[2] < e0[ie] .AND. e0[ie] < e[3]:
273+ elif e[2] < e0[ie] < e[3]:
271274 #
272275 v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e[0:4] - e0[ie])
273276 w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
--- a/src/libtetrabz/libtetrabz_occ.py
+++ b/src/libtetrabz/libtetrabz_occ.py
@@ -87,7 +87,7 @@ def occ(bvec=numpy.array([1.0, 0.0, 0.0]), eig=numpy.array([0.0])):
8787 w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
8888 for ib in range(nb):
8989 #
90- e = ei1[0:4, ib]
90+ e = ei1[0:4, ib].copy()
9191 indx = e.argsort(0)
9292 e.sort(0)
9393 #
--- a/src/libtetrabz/libtetrabz_polcmplx.py
+++ b/src/libtetrabz/libtetrabz_polcmplx.py
@@ -53,14 +53,15 @@ def polcmplx(bvec, eig1, eig2, e0):
5353 ei1 = wlsm.dot(eig1t)
5454 ej1 = wlsm.dot(eig2t)
5555 #
56- w1 = numpy.zeros([nb, nb, ne, 4], dtype=numpy.float_)
56+ w1 = numpy.zeros([nb, nb, ne, 4], dtype=numpy.complex_)
57+ w2t = numpy.zeros([nb, ne, 4], dtype=numpy.complex_)
5758 for ib in range(nb):
5859 #
59- e = ei1[0:4, ib]
60+ e = ei1[0:4, ib].copy()
6061 indx = e.argsort(0)
6162 e.sort(0)
6263 #
63- if e[0] <= 0.0 < e(2):
64+ if e[0] <= 0.0 < e[1]:
6465 #
6566 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
6667 #
@@ -68,9 +69,10 @@ def polcmplx(bvec, eig1, eig2, e0):
6869 ei2 = tsmall.dot(ei1[indx[0:4], ib])
6970 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7071 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
71- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
72+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
73+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
7274 #
73- elif e(2) <= 0.0 < e[2]:
75+ elif e[1] <= 0.0 < e[2]:
7476 #
7577 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
7678 #
@@ -78,7 +80,8 @@ def polcmplx(bvec, eig1, eig2, e0):
7880 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7981 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8082 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
81- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
83+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
84+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
8285 #
8386 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8487 #
@@ -86,7 +89,8 @@ def polcmplx(bvec, eig1, eig2, e0):
8689 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8790 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8891 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
89- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
92+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
93+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
9094 #
9195 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9296 #
@@ -94,7 +98,8 @@ def polcmplx(bvec, eig1, eig2, e0):
9498 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9599 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
96100 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
97- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
101+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
102+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
98103 #
99104 elif e[2] <= 0.0 < e[3]:
100105 #
@@ -104,7 +109,8 @@ def polcmplx(bvec, eig1, eig2, e0):
104109 ei2 = tsmall.dot(ei1[indx[0:4], ib])
105110 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
106111 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
107- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
112+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
113+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
108114 #
109115 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
110116 #
@@ -112,7 +118,8 @@ def polcmplx(bvec, eig1, eig2, e0):
112118 ei2 = tsmall.dot(ei1[indx[0:4], ib])
113119 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
114120 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
115- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
121+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
122+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
116123 #
117124 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
118125 #
@@ -120,7 +127,8 @@ def polcmplx(bvec, eig1, eig2, e0):
120127 ei2 = tsmall.dot(ei1[indx[0:4], ib])
121128 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
122129 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
123- w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
130+ w2t[0:nb, 0:ne, indx[0:4]] = w2.dot(tsmall)
131+ w1[ib, 0:nb, 0:ne, 0:4] += v * w2t
124132 #
125133 elif e[3] <= 0.0:
126134 #
@@ -151,6 +159,7 @@ def libtetrabz_polcmplx2(e0, ei1, ej1):
151159 nb = ej1.shape[1]
152160 thr = 1.0e-8
153161 w1 = numpy.zeros([nb, ne, 4], dtype=numpy.complex_)
162+ w2t = numpy.zeros([ne, 4], dtype=numpy.complex_)
154163
155164 for ib in range(nb):
156165 #
@@ -158,77 +167,69 @@ def libtetrabz_polcmplx2(e0, ei1, ej1):
158167 indx = e.argsort(0)
159168 e.sort(0)
160169 #
161- if e[0] <= 0.0 < e(2) or e[0] < 0.0 <= e[1]:
170+ if e[0] <= 0.0 < e[1] or e[0] < 0.0 <= e[1]:
162171 #
163172 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
164173 #
165174 if v > thr:
166- #
167175 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
168176 w2 = libtetrabz_polcmplx3(e0, de)
169- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
170- #
177+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
178+ w1[ib, 0:ne, 0:4] += v * w2t
171179 #
172180 elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
173181 #
174182 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
175183 #
176184 if v > thr:
177- #
178185 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
179186 w2 = libtetrabz_polcmplx3(e0, de)
180- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
181- #
187+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
188+ w1[ib, 0:ne, 0:4] += v * w2t
182189 #
183190 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
184191 #
185192 if v > thr:
186- #
187193 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
188194 w2 = libtetrabz_polcmplx3(e0, de)
189- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
190- #
195+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
196+ w1[ib, 0:ne, 0:4] += v * w2t
191197 #
192198 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
193199 #
194200 if v > thr:
195- #
196201 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
197202 w2 = libtetrabz_polcmplx3(e0, de)
198- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
199- #
203+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
204+ w1[ib, 0:ne, 0:4] += v * w2t
200205 #
201206 elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
202207 #
203208 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
204209 #
205210 if v > thr:
206- #
207211 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
208212 w2 = libtetrabz_polcmplx3(e0, de)
209- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
210- #
213+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
214+ w1[ib, 0:ne, 0:4] += v * w2t
211215 #
212216 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
213217 #
214218 if v > thr:
215- #
216219 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
217220 w2 = libtetrabz_polcmplx3(e0, de)
218- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
219- #
221+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
222+ w1[ib, 0:ne, 0:4] += v * w2t
220223 #
221224 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
222225 #
223226 if v > thr:
224- #
225227 de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
226228 w2 = libtetrabz_polcmplx3(e0, de)
227- w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
228- #
229+ w2t[0:ne, indx[0:4]] = w2.dot(tsmall)
230+ w1[ib, 0:ne, 0:4] += v * w2t
229231 #
230232 elif e[3] <= 0.0:
231- #
232233 de = ej1[0:4, ib] - ei1[0:4]
233234 w2 = libtetrabz_polcmplx3(e0, de)
234235 w1[ib, 0:ne, 0:4] += w2
@@ -244,7 +245,7 @@ def libtetrabz_polcmplx3(e0, de):
244245 :return:
245246 """
246247 ne = e0.shape[0]
247- e = de[0:4]
248+ e = de[0:4].copy()
248249 indx = e.argsort(0)
249250 e.sort(0)
250251 w1 = numpy.zeros([ne, 4], dtype=numpy.complex_)
@@ -252,7 +253,7 @@ def libtetrabz_polcmplx3(e0, de):
252253 #
253254 for ie in range(ne):
254255 #
255- # I don't know which one is better.
256+ # I am no sure which one is better.
256257 # The former is more stable.
257258 # The latter is more accurate ?
258259 #
@@ -262,11 +263,11 @@ def libtetrabz_polcmplx3(e0, de):
262263 #
263264 x = (e[0:4] + e0[ie].real) / e0[ie].imag
264265 # thr = maxval(de(1:4)) * 1d-3
265- thr = max(1.0e-3, x.max() * 1.0e-2)
266+ thr = max(1.0e-3, numpy.abs(x).max() * 1.0e-2)
266267 #
267- if abs(x[3] - x(3)) < thr:
268- if abs(x[3] - x(2)) < thr:
269- if abs(x[3] - x(1)) < thr:
268+ if abs(x[3] - x[2]) < thr:
269+ if abs(x[3] - x[1]) < thr:
270+ if abs(x[3] - x[0]) < thr:
270271 #
271272 # e[3] = e[2] = e[1] = e[0]
272273 #
@@ -280,10 +281,10 @@ def libtetrabz_polcmplx3(e0, de):
280281 #
281282 # e[3] = e[2] = e[1]
282283 #
283- w2[0:2, 3] = libtetrabz_polcmplx_1211(x[3], x(1))
284+ w2[0:2, 3] = libtetrabz_polcmplx_1211(x[3], x[0])
284285 w2[0:2, 2] = w2[0:2, 3]
285286 w2[0:2, 1] = w2[0:2, 3]
286- w2[0:2, 0] = libtetrabz_polcmplx_1222(x[0], x(4))
287+ w2[0:2, 0] = libtetrabz_polcmplx_1222(x[0], x[3])
287288 #
288289 # IF(ANY(w2(1:2,1:4) < 0.0)):
289290 # WRITE(*,*) ie
@@ -291,13 +292,13 @@ def libtetrabz_polcmplx3(e0, de):
291292 # WRITE(*,'(2e15.5)') w2(1:2,1:4)
292293 # STOP "weighting 4=3=2"
293294 #
294- elif abs(x[1] - x(1)) < thr:
295+ elif abs(x[1] - x[0]) < thr:
295296 #
296297 # e[3] = e[2], e[1] = e[0]
297298 #
298- w2[0:2, 3] = libtetrabz_polcmplx_1221(x[3], x(2))
299+ w2[0:2, 3] = libtetrabz_polcmplx_1221(x[3], x[1])
299300 w2[0:2, 2] = w2[0:2, 3]
300- w2[0:2, 1] = libtetrabz_polcmplx_1221(x[1], x(4))
301+ w2[0:2, 1] = libtetrabz_polcmplx_1221(x[1], x[3])
301302 w2[0:2, 0] = w2[0:2, 1]
302303 #
303304 # IF(ANY(w2(1:2,1:4) < 0.0)):
@@ -310,10 +311,10 @@ def libtetrabz_polcmplx3(e0, de):
310311 #
311312 # e[3] = e[2]
312313 #
313- w2[0:2, 3] = libtetrabz_polcmplx_1231(x[3], x[0], x(2))
314+ w2[0:2, 3] = libtetrabz_polcmplx_1231(x[3], x[0], x[1])
314315 w2[0:2, 2] = w2[0:2, 3]
315- w2[0:2, 1] = libtetrabz_polcmplx_1233(x[1], x[0], x(4))
316- w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[1], x(4))
316+ w2[0:2, 1] = libtetrabz_polcmplx_1233(x[1], x[0], x[3])
317+ w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[1], x[3])
317318 #
318319 # IF(ANY(w2(1:2,1:4) < 0.0)):
319320 # WRITE(*,*) ie
@@ -321,13 +322,13 @@ def libtetrabz_polcmplx3(e0, de):
321322 # WRITE(*,'(2e15.5)') w2(1:2,1:4)
322323 # STOP "weighting 4=3"
323324 #
324- elif abs(x[2] - x(2)) < thr:
325- if abs(x[2] - x(1)) < thr:
325+ elif abs(x[2] - x[1]) < thr:
326+ if abs(x[2] - x[0]) < thr:
326327 #
327328 # e[2] = e[1] = e[0]
328329 #
329- w2[0:2, 3] = libtetrabz_polcmplx_1222(x[3], x(3))
330- w2[0:2, 2] = libtetrabz_polcmplx_1211(x[2], x(4))
330+ w2[0:2, 3] = libtetrabz_polcmplx_1222(x[3], x[2])
331+ w2[0:2, 2] = libtetrabz_polcmplx_1211(x[2], x[3])
331332 w2[0:2, 1] = w2[0:2, 2]
332333 w2[0:2, 0] = w2[0:2, 2]
333334 #
@@ -341,10 +342,10 @@ def libtetrabz_polcmplx3(e0, de):
341342 #
342343 # e[2] = e[1]
343344 #
344- w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[0], x(3))
345- w2[0:2, 2] = libtetrabz_polcmplx_1231(x[2], x[0], x(4))
345+ w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[0], x[2])
346+ w2[0:2, 2] = libtetrabz_polcmplx_1231(x[2], x[0], x[3])
346347 w2[0:2, 1] = w2[0:2, 2]
347- w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[3], x(3))
348+ w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[3], x[2])
348349 #
349350 # IF(ANY(w2(1:2,1:4) < 0.0)):
350351 # WRITE(*,*) ie
@@ -352,13 +353,13 @@ def libtetrabz_polcmplx3(e0, de):
352353 # WRITE(*,'(2e15.5)') w2(1:2,1:4)
353354 # STOP "weighting 3=2"
354355 #
355- elif abs(x[1] - x(1)) < thr:
356+ elif abs(x[1] - x[0]) < thr:
356357 #
357358 # e[1] = e[0]
358359 #
359- w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[2], x(2))
360- w2[0:2, 2] = libtetrabz_polcmplx_1233(x[2], x[3], x(2))
361- w2[0:2, 1] = libtetrabz_polcmplx_1231(x[1], x[2], x(4))
360+ w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[2], x[1])
361+ w2[0:2, 2] = libtetrabz_polcmplx_1233(x[2], x[3], x[1])
362+ w2[0:2, 1] = libtetrabz_polcmplx_1231(x[1], x[2], x[3])
362363 w2[0:2, 0] = w2[0:2, 1]
363364 #
364365 # IF(ANY(w2(1:2,1:4) < 0.0)):
@@ -371,10 +372,10 @@ def libtetrabz_polcmplx3(e0, de):
371372 #
372373 # Different each other.
373374 #
374- w2[0:2, 3] = libtetrabz_polcmplx_1234(x[3], x[0], x[1], x(3))
375- w2[0:2, 2] = libtetrabz_polcmplx_1234(x[2], x[0], x[1], x(4))
376- w2[0:2, 1] = libtetrabz_polcmplx_1234(x[1], x[0], x[2], x(4))
377- w2[0:2, 0] = libtetrabz_polcmplx_1234(x[0], x[1], x[2], x(4))
375+ w2[0:2, 3] = libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2])
376+ w2[0:2, 2] = libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3])
377+ w2[0:2, 1] = libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3])
378+ w2[0:2, 0] = libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3])
378379 #
379380 # IF(ANY(w2(1:2,1:4) < 0.0)):
380381 # WRITE(*,*) ie
@@ -382,9 +383,9 @@ def libtetrabz_polcmplx3(e0, de):
382383 # WRITE(*,'(2e15.5)') w2(1:2,1:4)
383384 # STOP "weighting"
384385 #
385- #
386- w1[ie, indx[0:4]] = complex(w2[0, 0:4] / e0[ie].imag, w2[1, 0:4] / (- e0[ie].imag))
387- #
386+ #
387+ w1[ie, indx[0:4]] = w2[0, 0:4] / e0[ie].imag + 1.0j * w2[1, 0:4] / (- e0[ie].imag)
388+ #
388389 return w1
389390 #
390391 # Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
--- a/src/libtetrabz/libtetrabz_polstat.py
+++ b/src/libtetrabz/libtetrabz_polstat.py
@@ -53,9 +53,10 @@ def polstat(bvec, eig1, eig2):
5353 ej1 = wlsm.dot(eig2t)
5454 #
5555 w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
56+ w2t = numpy.zeros([nb, 4], dtype=numpy.float_)
5657 for ib in range(nb):
5758 #
58- e = ei1[0:4, ib]
59+ e = ei1[0:4, ib].copy()
5960 indx = e.argsort(0)
6061 e.sort(0)
6162 #
@@ -67,7 +68,8 @@ def polstat(bvec, eig1, eig2):
6768 ei2 = tsmall.dot(ei1[indx[0:4], ib])
6869 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
6970 w2 = libtetrabz_polstat2(ei2, ej2)
70- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
71+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
72+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
7173 #
7274 elif e[1] <= 0.0 < e[2]:
7375 #
@@ -77,7 +79,8 @@ def polstat(bvec, eig1, eig2):
7779 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7880 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7981 w2 = libtetrabz_polstat2(ei2, ej2)
80- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
82+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
83+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
8184 #
8285 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8386 #
@@ -85,7 +88,8 @@ def polstat(bvec, eig1, eig2):
8588 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8689 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8790 w2 = libtetrabz_polstat2(ei2, ej2)
88- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
91+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
92+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
8993 #
9094 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9195 #
@@ -93,7 +97,8 @@ def polstat(bvec, eig1, eig2):
9397 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9498 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9599 w2 = libtetrabz_polstat2(ei2, ej2)
96- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
100+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
101+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
97102 #
98103 elif e[2] <= 0.0 < e[3]:
99104 #
@@ -103,7 +108,8 @@ def polstat(bvec, eig1, eig2):
103108 ei2 = tsmall.dot(ei1[indx[0:4], ib])
104109 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
105110 w2 = libtetrabz_polstat2(ei2, ej2)
106- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
111+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
112+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
107113 #
108114 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
109115 #
@@ -111,7 +117,8 @@ def polstat(bvec, eig1, eig2):
111117 ei2 = tsmall.dot(ei1[indx[0:4], ib])
112118 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
113119 w2 = libtetrabz_polstat2(ei2, ej2)
114- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
120+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
121+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
115122 #
116123 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
117124 #
@@ -119,7 +126,8 @@ def polstat(bvec, eig1, eig2):
119126 ei2 = tsmall.dot(ei1[indx[0:4], ib])
120127 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
121128 w2 = libtetrabz_polstat2(ei2, ej2)
122- w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
129+ w2t[0:nb, indx[0:4]] = w2.dot(tsmall)
130+ w1[ib, 0:nb, 0:4] += v * w2t[0:nb, 0:4]
123131 #
124132 elif e[3] <= 0.0:
125133 ei2 = ei1[0:4, ib]
@@ -134,11 +142,12 @@ def polstat(bvec, eig1, eig2):
134142 wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb] += w1.dot(wlsm[:, ii])
135143 #
136144 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
145+ return wght
137146
138147
139148 def libtetrabz_polstat2(ei1, ej1):
140149 """
141- Tetrahedra method for theta( - E2)
150+ Tetrahedron method for theta( - E2)
142151 :param ei1:
143152 :param ej1:
144153 :return:
@@ -231,7 +240,7 @@ def libtetrabz_polstat3(de):
231240 :return:
232241 """
233242 #
234- e = de[0:4]
243+ e = de[0:4].copy()
235244 indx = e.argsort(0)
236245 e.sort(0)
237246 #
@@ -274,7 +283,7 @@ def libtetrabz_polstat3(de):
274283 print(w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]])
275284 raise ValueError("weighting 4=3=2")
276285 #
277- elif abs(e[1] - e(1)) < thr:
286+ elif abs(e[1] - e[0]) < thr:
278287 #
279288 # e[3] = e[2], e[1] = e[0]
280289 #
@@ -302,8 +311,8 @@ def libtetrabz_polstat3(de):
302311 print(w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]])
303312 raise ValueError("weighting 4=3")
304313 #
305- elif abs(e[2] - e(2)) < thr:
306- if abs(e[2] - e(1)) < thr:
314+ elif abs(e[2] - e[1]) < thr:
315+ if abs(e[2] - e[0]) < thr:
307316 #
308317 # e[2] = e[1] = e[0]
309318 #
@@ -331,7 +340,7 @@ def libtetrabz_polstat3(de):
331340 print(w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]])
332341 raise ValueError("weighting 3=2")
333342 #
334- elif abs(e[1] - e(1)) < thr:
343+ elif abs(e[1] - e[0]) < thr:
335344 #
336345 # e[1] = e[0]
337346 #
--- /dev/null
+++ b/src/lindhard.py
@@ -0,0 +1,57 @@
1+import numpy
2+import libtetrabz
3+import matplotlib.pyplot as plt
4+import matplotlib
5+
6+bvec = numpy.array([[3.0, 0.0, 0.0],
7+ [0.0, 3.0, 0.0],
8+ [0.0, 0.0, 3.0]])
9+VBZ = abs(numpy.linalg.det(bvec))
10+
11+ng0 = 10
12+nb = 1
13+ng = numpy.array([ng0, ng0, ng0])
14+eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
15+
16+ef = 0.5
17+qmax = 4.0
18+
19+eig1 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
20+eig2 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
21+qx = numpy.arange(0.0, qmax, 0.1)
22+chi = numpy.empty(qx.shape, dtype=numpy.float_)
23+for iq in range(qx.shape[0]):
24+ print(iq, "in", qx.shape[0])
25+ #
26+ qvec = numpy.array([qx[iq], 0.0, 0.0])
27+ #
28+ for i0 in range(ng[0]):
29+ for i1 in range(ng[1]):
30+ for i2 in range(ng[2]):
31+ #
32+ ikvec = numpy.array([i0, i1, i2])
33+ for ii in range(3):
34+ if ikvec[ii] * 2 >= ng[ii]:
35+ ikvec[ii] += - ng[ii]
36+ kvec = ikvec[0:3] / ng[0:3]
37+ kvec[0:3] = kvec.dot(bvec)
38+ #
39+ eig1[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec) - ef
40+ #
41+ kvec[0:3] += qvec[0:3]
42+ eig2[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec) - ef
43+ #
44+ #
45+ if iq == 0:
46+ e0 = numpy.array([0.0])
47+ wght = libtetrabz.dos(bvec, eig1, e0)
48+ chi[iq] = wght.sum() * VBZ / (4.0*numpy.pi)
49+ else:
50+ wght = libtetrabz.polstat(bvec, eig1, eig2)
51+ chi[iq] = wght.sum() * VBZ / (4.0 * numpy.pi) * 2.0
52+ #
53+matplotlib.use('TkAgg')
54+chi0 = 0.5+0.5/qx*(1-0.25*qx**2)*numpy.log(numpy.abs((qx+2)/(qx-2)))
55+plt.plot(qx, chi, label="Calc")
56+plt.plot(qx, chi0, label="Exact")
57+plt.show()
--- /dev/null
+++ b/src/test.py
@@ -0,0 +1,359 @@
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 numpy
25+import libtetrabz
26+
27+
28+def test_occ(nb, bvec, vbz, eig, mat):
29+ """
30+
31+ :return:
32+ """
33+ #
34+ wght = libtetrabz.occ(bvec, eig - 0.5)
35+ #
36+ val = numpy.empty(2, dtype=numpy.float_)
37+ for ib in range(nb):
38+ val[ib] = numpy.sum(wght[:, :, :, ib] * mat[:, :, :])
39+ print("# libtetrabz_occ")
40+ print(' %15.5e %15.5e' % (4.0 * numpy.pi / 5.0, val[0] * vbz))
41+ print(' %15.5e %15.5e' % (numpy.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))
42+ print("")
43+
44+
45+def test_fermieng(nb, bvec, vbz, eig, mat):
46+ """
47+
48+ :return:
49+ """
50+ #
51+ nelec = (4.0 * numpy.pi / 3.0 + math.sqrt(2.0) * numpy.pi / 3.0) / vbz
52+
53+ ef, wght, iterations = libtetrabz.fermieng(bvec, eig, nelec)
54+ #
55+ val = numpy.empty(2, dtype=numpy.float_)
56+ for ib in range(nb):
57+ val[ib] = numpy.sum(wght[:, :, :, ib] * mat[:, :, :])
58+ #
59+ print("# libtetrabz_fermieng")
60+ print(" %15.5e %15.5e" % (0.5, ef))
61+ print(" %15.5e %15.5e" % (4.0 * numpy.pi / 5.0, val[0] * vbz))
62+ print(" %15.5e %15.5e" % (numpy.pi / (5.0 * math.sqrt(2.0)), val[1] * vbz))
63+ print("")
64+
65+
66+def test_dos(nb, bvec, vbz, eig, mat):
67+ """
68+
69+ :return:
70+ """
71+ ne = 5
72+ #
73+ e0 = numpy.empty(ne, dtype=numpy.float_)
74+ val0 = numpy.empty([nb, ne], dtype=numpy.float_)
75+ val = numpy.empty([nb, ne], dtype=numpy.float_)
76+ for ie in range(ne):
77+ e0[ie] = 0.2 * (ie + 1)
78+ val0[0, ie] = 4.0 * numpy.pi * e0[ie]**3
79+ if e0[ie] > 1.0 / math.sqrt(2.0):
80+ val0[1, ie] = math.sqrt(2.0) * numpy.pi * math.sqrt(-1.0 + 2.0 * e0[ie]**2)**3
81+ else:
82+ val0[1, ie] = 0.0
83+ e0[0:ne] = e0[0:ne]**2 * 0.5
84+ #
85+ wght = libtetrabz.dos(bvec, eig, e0)
86+ #
87+ for ib in range(nb):
88+ for ie in range(ne):
89+ val[ib, ie] = numpy.sum(wght[:, :, :, ib, ie] * mat[:, :, :])
90+ #
91+ print("# libtetrabz_dos")
92+ for ib in range(nb):
93+ for ie in range(ne):
94+ print(" %15.5e %15.5e" % (val0[ib, ie], val[ib, ie] * vbz))
95+ print("")
96+
97+
98+def test_intdos(nb, bvec, vbz, eig, mat):
99+ """
100+
101+ :return:
102+ """
103+ ne = 5
104+ e0 = numpy.empty(ne, dtype=numpy.float_)
105+ val0 = numpy.empty([nb, ne], dtype=numpy.float_)
106+ val = numpy.empty([nb, ne], dtype=numpy.float_)
107+ #
108+ for ie in range(ne):
109+ e0[ie] = 0.2 * (ie + 1)
110+ val0[0, ie] = 4.0 * numpy.pi * e0[ie]**5 / 5.0
111+ if e0[ie] > 1.0 / math.sqrt(2.0):
112+ val0[1, ie] = numpy.pi * math.sqrt(-1.0 + 2.0 * e0[ie]**2)**5 / (5.0 * math.sqrt(2.0))
113+ else:
114+ val0[1, ie] = 0.0
115+ e0[0:ne] = e0[0:ne]**2 * 0.5
116+ #
117+ wght = libtetrabz.intdos(bvec, eig, e0)
118+ #
119+ for ib in range(nb):
120+ for ie in range(ne):
121+ val[ib, ie] = numpy.sum(wght[:, :, :, ib, ie] * mat[:, :, :])
122+ #
123+ print("# libtetrabz_intdos")
124+ for ib in range(nb):
125+ for ie in range(ne):
126+ print(" %15.5e %15.5e" % (val0[ib, ie], val[ib, ie] * vbz))
127+ print("")
128+
129+
130+def test_dblstep(nb, bvec, vbz, eig1, eig2, mat):
131+ """
132+
133+ :param bvec:
134+ :param vbz:
135+ :param eig1:
136+ :param eig2:
137+ :param mat:
138+ :return:
139+ """
140+ #
141+ wght = libtetrabz.dblstep(bvec, eig1 - 0.5, eig2 - 0.5)
142+ #
143+ val = numpy.empty([nb, nb], dtype=numpy.float_)
144+ for ib in range(nb):
145+ for jb in range(nb):
146+ val[ib, jb] = numpy.sum(wght[:, :, :, ib, jb] * mat[:, :, :])
147+ #
148+ print("# libtetrabz_dblstep")
149+ print(" %15.5e %15.5e" % (49.0 * numpy.pi / 320.0, val[0, 0] * vbz))
150+ print(" %15.5e %15.5e" % (0.0, val[0, 1] * vbz))
151+ print(" %15.5e %15.5e" % (numpy.pi * (512.0 * math.sqrt(2.0) - 319.0) / 10240.0, val[1, 0] * vbz))
152+ print(" %15.5e %15.5e" % (0.0, val[1, 1] * vbz))
153+ print("")
154+
155+
156+def test_dbldelta(nb, bvec, vbz, eig1, eig2, mat):
157+ """
158+
159+ :param bvec:
160+ :param vbz:
161+ :param eig1:
162+ :param eig2:
163+ :param mat:
164+ :return:
165+ """
166+ #
167+ wght = libtetrabz.dbldelta(bvec, eig1 - 0.5, eig2 - 0.5)
168+ #
169+ val = numpy.empty([nb, nb], dtype=numpy.float_)
170+ for ib in range(nb):
171+ for jb in range(nb):
172+ val[ib, jb] = numpy.sum(wght[:, :, :, ib, jb] * mat[:, :, :])
173+ print("# libtetrabz_dbldelta")
174+ print(" %15.5e %15.5e" % (2.0 * numpy.pi, val[0, 0] * vbz))
175+ print(" %15.5e %15.5e" % (0.0, val[0, 1] * vbz))
176+ print(" %15.5e %15.5e" % (numpy.pi, val[1, 0] * vbz))
177+ print(" %15.5e %15.5e" % (0.0, val[1, 1] * vbz))
178+ print("")
179+
180+
181+def test_polstat(nb, bvec, vbz, eig1, eig2, mat):
182+ """
183+
184+ :param bvec:
185+ :param vbz:
186+ :param eig1:
187+ :param eig2:
188+ :param mat:
189+ :return:
190+ """
191+ wght = libtetrabz.polstat(bvec, eig1 - 0.5, eig2 - 0.5)
192+ #
193+ val = numpy.empty([nb, nb], dtype=numpy.float_)
194+ for ib in range(nb):
195+ for jb in range(nb):
196+ val[ib, jb] = numpy.sum(wght[:, :, :, ib, jb] * mat[:, :, :])
197+ #
198+ print("# libtetrabz_polstat")
199+ print(" %15.5e %15.5e" % (numpy.pi * (68.0 + 45.0 * math.log(3.0)) / 96.0, val[0, 0] * vbz))
200+ print(" %15.5e %15.5e" % (numpy.pi * 8.0 / 5.0, val[0, 1] * vbz))
201+ print(" %15.5e %15.5e" % (numpy.pi * (228.0 + 22.0 * math.sqrt(2.0) - 96.0 * math.log(2.0)
202+ + 192.0 * math.log(4.0 + math.sqrt(2.0))
203+ - 3.0 * math.log(1.0 + 2.0 * math.sqrt(2.0))) / 1536.0,
204+ val[1, 0] * vbz))
205+ print(" %15.5e %15.5e" % (numpy.pi * math.sqrt(8.0) / 5.0, val[1, 1] * vbz))
206+ print("")
207+
208+
209+def test_fermigr(nb, bvec, vbz, eig1, eig2, mat):
210+ """
211+
212+ :param bvec:
213+ :param vbz:
214+ :param eig1:
215+ :param eig2:
216+ :param mat:
217+ :return:
218+ """
219+ ne = 3
220+ #
221+ e0 = numpy.empty(ne, dtype=numpy.float_)
222+ val0 = numpy.empty([nb, nb, ne], dtype=numpy.float_)
223+ val = numpy.empty([nb, nb, ne], dtype=numpy.float_)
224+ e0[0] = 1.0 / 3.0
225+ e0[1] = 2.0 / 3.0
226+ e0[2] = 1.0
227+ val0[0, 0, 0] = 4.0 * numpy.pi / 9.0
228+ val0[0, 0, 1] = 1295.0 * numpy.pi / 2592.0
229+ val0[0, 0, 2] = 15.0 * numpy.pi / 32.0
230+ val0[1, 0, 0] = 5183.0 * numpy.pi / 41472.0
231+ val0[1, 0, 1] = 4559.0 * numpy.pi / 41472.0
232+ val0[1, 0, 2] = 0.0
233+ val0[0:nb, 1, 0:3] = 0.0
234+ #
235+ wght = libtetrabz.fermigr(bvec, eig1 - 0.5, eig2 - 0.5, e0)
236+ #
237+ for ib in range(nb):
238+ for jb in range(nb):
239+ for ie in range(ne):
240+ val[ib, jb, ie] = numpy.sum(wght[:, :, :, ib, jb, ie] * mat[:, :, :])
241+ #
242+ print("# libtetrabz_fermigr")
243+ for ib in range(nb):
244+ for jb in range(nb):
245+ for ie in range(ne):
246+ print(" %15.5e %15.5e" % (val0[ib, jb, ie], val[ib, jb, ie] * vbz))
247+ print("")
248+
249+
250+def test_polcmplx(nb, bvec, vbz, eig1, eig2, mat):
251+ """
252+
253+ :param bvec:
254+ :param vbz:
255+ :param eig1:
256+ :param eig2:
257+ :param mat:
258+ :return:
259+ """
260+ ne = 3
261+ e0 = numpy.empty(ne, dtype=numpy.complex_)
262+ val0 = numpy.empty([nb, nb, ne], dtype=numpy.complex_)
263+ val = numpy.empty([nb, nb, ne], dtype=numpy.complex_)
264+ e0[0] = -2.0 + 1.0j
265+ e0[1] = 0.0 + 2.0j
266+ e0[2] = 1.0 - 0.5j
267+ val0[0, 0, 0] = -0.838243341280338 - 0.734201894333234j
268+ val0[0, 0, 1] = 0.270393588876530 - 0.771908416949610j
269+ val0[0, 0, 2] = 0.970996830573510 + 0.302792326476720j
270+ val0[1, 0, 0] = -0.130765724778920 - 0.087431218706638j
271+ val0[1, 0, 1] = 0.030121954547245 - 0.135354254293510j
272+ val0[1, 0, 2] = 0.178882244951203 + 0.064232167683425j
273+ val0[0, 1, 0:3] = (8.0 * numpy.pi) / (5.0 * (1.0 + 2.0 * e0[0:3]))
274+ val0[1, 1, 0:3] = (math.sqrt(8.0) * numpy.pi) / (5.0 * (1.0 + 4.0 * e0[0:3]))
275+ #
276+ wght = libtetrabz.polcmplx(bvec, eig1 - 0.5, eig2 - 0.5, e0)
277+ #
278+ for ib in range(nb):
279+ for jb in range(nb):
280+ for ie in range(ne):
281+ val[ib, jb, ie] = numpy.sum(wght[:, :, :, ib, jb, ie] * mat[:, :, :])
282+ #
283+ print("# libtetrabz_polcmplx")
284+ for ib in range(nb):
285+ for jb in range(nb):
286+ for ie in range(ne):
287+ print(" %15.5e %15.5e" % (val0[ib, jb, ie].real, val[ib, jb, ie].real * vbz))
288+ print(" %15.5e %15.5e" % (val0[ib, jb, ie].imag, val[ib, jb, ie].imag * vbz))
289+ print("")
290+
291+
292+def test():
293+ ng0 = 8
294+ ng = numpy.array([ng0, ng0, ng0])
295+ nb = 2
296+ bvec = numpy.array([[3.0, 0.0, 0.0],
297+ [0.0, 3.0, 0.0],
298+ [0.0, 0.0, 3.0]])
299+ vbz = abs(numpy.linalg.det(bvec))
300+
301+ #
302+ eig1 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
303+ eig2 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
304+ mat = numpy.empty([ng[0], ng[1], ng[2]], dtype=numpy.float_)
305+ #
306+ kvec = numpy.empty(3, dtype=numpy.float_)
307+ for i0 in range(ng[1]):
308+ for i1 in range(ng[1]):
309+ for i2 in range(ng[2]):
310+ #
311+ kvec[0:3] = numpy.array([i0, i1, i2]) / ng[0:3]
312+ for ii in range(3):
313+ if kvec[ii] > 0.5:
314+ kvec[ii] += -1
315+ kvec[0:3] = kvec.dot(bvec)
316+ #
317+ eig1[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec)
318+ eig1[i0, i1, i2, 1] = eig1[i0, i1, i2, 0] + 0.25
319+ #
320+ kvec[0] += 1.0
321+ eig2[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec)
322+ eig2[i0, i1, i2, 1] = eig1[i0, i1, i2, 0] + 0.5
323+ #
324+ #
325+ for i0 in range(ng[1]):
326+ for i1 in range(ng[1]):
327+ for i2 in range(ng[2]):
328+ #
329+ kvec[0:3] = numpy.array([i0, i1, i2]) / ng[0:3]
330+ for ii in range(3):
331+ if kvec[ii] > 0.5:
332+ kvec[ii] += -1
333+ kvec[0:3] = kvec.dot(bvec)
334+ #
335+ mat[i0, i1, i2] = kvec.dot(kvec)
336+ #
337+ print("# Ideal Result")
338+ #
339+ test_occ(nb, bvec, vbz, eig1, mat)
340+ #
341+ test_fermieng(nb, bvec, vbz, eig1, mat)
342+ #
343+ test_dos(nb, bvec, vbz, eig1, mat)
344+ #
345+ test_intdos(nb, bvec, vbz, eig1, mat)
346+ #
347+ test_dblstep(nb, bvec, vbz, eig1, eig2, mat)
348+ #
349+ test_dbldelta(nb, bvec, vbz, eig1, eig2, mat)
350+ #
351+ test_polstat(nb, bvec, vbz, eig1, eig2, mat)
352+ #
353+ test_fermigr(nb, bvec, vbz, eig1, eig2, mat)
354+ #
355+ test_polcmplx(nb, bvec, vbz, eig1, eig2, mat)
356+ #
357+
358+
359+test()
--- /dev/null
+++ b/src/tutorial.ipynb
@@ -0,0 +1,254 @@
1+{
2+ "cells": [
3+ {
4+ "cell_type": "code",
5+ "execution_count": 1,
6+ "id": "4cdce61d",
7+ "metadata": {},
8+ "outputs": [
9+ {
10+ "name": "stderr",
11+ "output_type": "stream",
12+ "text": [
13+ "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib\\site-packages\\numpy\\_distributor_init.py:32: UserWarning: loaded more than 1 DLL from .libs:\n",
14+ "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib\\site-packages\\numpy\\.libs\\libopenblas.WCDJNK7YVMPZQ2ME2ZZHJJRJ3JIKNDB7.gfortran-win_amd64.dll\n",
15+ "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib\\site-packages\\numpy\\.libs\\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll\n",
16+ " stacklevel=1)\n"
17+ ]
18+ }
19+ ],
20+ "source": [
21+ "import numpy\n",
22+ "import libtetrabz\n",
23+ "import matplotlib.pyplot as plt"
24+ ]
25+ },
26+ {
27+ "cell_type": "code",
28+ "execution_count": 2,
29+ "id": "43f60703",
30+ "metadata": {},
31+ "outputs": [],
32+ "source": [
33+ "bvec=numpy.array([[10.0,0.0,0.0],\n",
34+ " [0.0,10.0,0.0],\n",
35+ " [0.0,0.0,10.0]])"
36+ ]
37+ },
38+ {
39+ "cell_type": "code",
40+ "execution_count": 3,
41+ "id": "bb4b70ad",
42+ "metadata": {},
43+ "outputs": [],
44+ "source": [
45+ "ng0 = 10\n",
46+ "nb = 1\n",
47+ "ng = numpy.array([ng0, ng0, ng0])\n",
48+ "eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
49+ "for i0 in range(ng[0]):\n",
50+ " for i1 in range(ng[1]):\n",
51+ " for i2 in range(ng[2]):\n",
52+ " kvec = 2.0*numpy.pi*numpy.array([i0, i1, i2]) / ng[0:3]\n",
53+ " eig[i0, i1, i2, 0] = - numpy.cos(kvec).sum(0)"
54+ ]
55+ },
56+ {
57+ "cell_type": "code",
58+ "execution_count": 4,
59+ "id": "4ffb1f01",
60+ "metadata": {},
61+ "outputs": [],
62+ "source": [
63+ "e0=numpy.linspace(-3, 3, 100)"
64+ ]
65+ },
66+ {
67+ "cell_type": "code",
68+ "execution_count": 5,
69+ "id": "7bc36b9c",
70+ "metadata": {
71+ "scrolled": false
72+ },
73+ "outputs": [],
74+ "source": [
75+ "wght=libtetrabz.dos(bvec,eig,e0)"
76+ ]
77+ },
78+ {
79+ "cell_type": "code",
80+ "execution_count": 6,
81+ "id": "75396ee0",
82+ "metadata": {},
83+ "outputs": [],
84+ "source": [
85+ "dos=wght.sum(3).sum(2).sum(1).sum(0)"
86+ ]
87+ },
88+ {
89+ "cell_type": "code",
90+ "execution_count": 7,
91+ "id": "72caf309",
92+ "metadata": {},
93+ "outputs": [
94+ {
95+ "data": {
96+ "text/plain": [
97+ "[<matplotlib.lines.Line2D at 0x160bd035c08>]"
98+ ]
99+ },
100+ "execution_count": 7,
101+ "metadata": {},
102+ "output_type": "execute_result"
103+ },
104+ {
105+ "data": {
106+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsSElEQVR4nO3deXhU5d3/8fd3ZrKvQBIISSAEwr4n4ErBHWwLbggoSqsW0Wpd+rOly1Ottk8X2+rTalVafOqOiBu2WuqCCwpCAmEJmyEBkrAkLNnXmdy/PzL4BAxmgJmcWb6v68plZuacyecY8snJOfe5jxhjUEopFbxsVgdQSinlW1r0SikV5LTolVIqyGnRK6VUkNOiV0qpIOewOsCJkpKSTGZmptUxlFIqoOTn5x8yxiR39prfFX1mZiZ5eXlWx1BKqYAiIntO9ppHh25EZKqI7BCRIhFZ2MnrC0Rks4gUiMgqERne4bWfuNfbISKXnd4mKKWUOl1dFr2I2IHHgWnAcGBOxyJ3e9EYM8oYMxb4PfAn97rDgdnACGAq8Ff3+ymllOomnuzRTwSKjDHFxpgWYAkwo+MCxpiaDg9jgGOX284Alhhjmo0xJUCR+/2UUkp1E0+O0acBpR0elwFnnbiQiHwfuBcIBy7ssO6aE9ZN62Td+cB8gH79+nmSWymllIe8NrzSGPO4MWYg8GPg56e47iJjTK4xJjc5udOTxkoppU6TJ0VfDmR0eJzufu5klgBXnOa6SimlvMyTol8HZIvIABEJp/3k6vKOC4hIdoeH3wS+cH++HJgtIhEiMgDIBtaeeWyllFKe6vIYvTHGKSJ3ACsAO/C0MaZQRB4E8owxy4E7RORioBU4Csxzr1soIkuBrYAT+L4xxuWjbVF+ztVm2La/hvKqRvZXNRIXGcZV49MQEaujBaWqhhaW5ZcRG+EgNTGKzF7R9O8VY3UsZQHxt/noc3NzjV4wFZz+sGIHj60sOu65BZMHsnDaUIsSBa+6ZifX/20NG8uqj3v+9dvPZVy/HhalUr4kIvnGmNzOXvO7K2NV8Pq85DBD+8Tx8DVjSE2M5NH3dvLkR7tIiArjtikDvfI1qhpaiAq3E+EIrMs1qhtbiXDYiAw789xNrS7mP5vHln01PHVDDiP6xlNyqJ4bFq9lbckRLfoQpEWvuoXT1caW8hpmTchgVHoCAA9OH0lNo5Pf/Xs7cZEOrj+rn0eHcWqbWtl+oJaoMDvJcRHER4bxwfYKluaV8vEXlUSH2Zk8JJmLh/XmkuG9iYsM8/XmnZa9hxv4d+F+3t16kPw9R4kOd/DtMX25NjedYanxVNY2U1HbRITDzpA+cYTZuz6l1uJs464lG/hs12EemTWGy0b0ASC9RzTpPaLYdMIevgoNWvSqW+yqrKex1cVod8kD2GzCH68dQ21TKz9/Ywt//6SYaaNSOX9QEhW1TeyqqKf0aAMADpuNNmPYuq+GnRW1dHbEMTUhkgWTB1LV0Mr72w7y9uYDxEU4uOGc/tx0/gCSYiNOKbPT1UbxoXq27qthz+EGjPs6wHCHjcEpcYxIi6dPfOQpn2PYUl7NXz8s4p0tBzAGhqfG8/0LBrGvqok3NpTz0tq9X1knMszG6LREMpOicbYZWl2GCIeNrOQYBibHIsCKwoO8u/UANU1OHvj2cK4cl37ce4xOT2BTedUpZVXBQYtedYuNZVUAjE5PPO75MLuNJ+bm8PqGct7evJ9FHxfzxIe7ALDbhL6JkdhEcLoMxhgG94nj8lGpjEqPp9VlqKht5khdC+P6JXLeoCTstvbSbWsbyYbSKhavKuaJj3axeFUJ5wzsRWpCFKkJkQDsr25kX1UTNU2tGNN+ObfT1UZjq4vGFhdH6ltodrZ97XbFRTqIjwwjOtxOVLgdEUHc2XvFhNM3MYre8ZHUNbeyv7qJkkP1bNhbRVyEg9smD+S6s/qR3iP6y/d7YPpw3tl8gMq6ZpLjIkiJi6CmyUnB3irW7z3KxzsP4bALYXYbDS1OluWXfblufKSDS4b34YpxfZmU/dXrUUanJ/L25gMcqW+hZ0z4KX4HVSDTolfdYlNZe7llJX111EdkmJ05E/sxZ2I/jta3sKm8mrTEKPr1jCbccXrX9NlsQk7/HuT0z2FXZR1//6SEzeVVbCmv5lBdCwBJseGkJkSRGN1+aEdEcNiEqHA70WF2esSEMyw1juGpCWQlx+Bw/xKpb3Gx40ANhftq+OJgHfUtThpbXDS2umgzYIzB1WYoOVTPZ7sOU9fsxGETesdH0ichkh9NHcLcs/sT38khpbjIMK6dkPGV56eP6dvpdtY2tVJyqJ76Zhc5/Xt87f+vY39NbSqrYsqQlFP7H6oCmha96habyqoZmZaAzfb1hzl6xIQzebB3r44emBzLb64a9eXjZqcLYzjtE5+xEQ5y+vckp39Pj5avb3YSGWb/8q8Nb4qLDPvKX0knMyotAZH274UWfWjRO0wpn2t2uti2v4bRGQldL9wNIhx2r4xu8VRMhMMnJX+q4iLDyEqK0ROyIUiLXvncjgO1tLoMo9MSrY4S8kanJ7LJfb5EhQ4teuVzxy7a6TjiRlljdHoCFbXNHKhusjqK6kZa9MrnNpVW0TMmnPQeUVZHCXnHjudv1L36kKJFr3xuc3k1o9MTdE4bPzCibzx2m7BZj9OHFC165VMNLU52HqxldJoetvEHkWF2BveO0z36EKNFr3yqcF8NbearF0op64xJT2BzeTX+NqGh8h0teuVTG0urAPxmaKVq/6Vb1dDK3iMNVkdR3USLXvlU4b4aesdHkBIXaXUU5XZs9NOW8hqLk6juokWvfGpXZR3ZKXFWx1AdZCW3T0Oxq7LO4iSqu2jRK58xxlBcWf9lsSj/EB3uoG9CJMVa9CFDi175TGVtM3XNzk4nMlPWykqOpfhQvdUxVDfRolc+c6xIspJjLU6iTpSVHENJZb2OvAkRWvTKZ4orjxW97tH7m6ykGGqbnVTWNVsdRXUDLXrlM8WVdUSG2eiboFMf+Jtjf2Ud+2WsgpsWvfKZ4kP1ZPaK6XIOetX9jv2VpUUfGrTolc8UV9YxUI/P+6W+CVFEhtl05E2I0KJXPtHibKP0aKMen/dTNpuQ2StGR96ECC165RN7j9TjajMM0KGVfisrOUb36EOER0UvIlNFZIeIFInIwk5ev1dEtorIJhF5X0T6d3jNJSIF7o/l3gyv/NeuSh1a6e+ykmIpPdpIi7PN6ijKx7osehGxA48D04DhwBwRGX7CYhuAXGPMaGAZ8PsOrzUaY8a6P6Z7Kbfyczq00v9lJcfgajPsPaKHb4KdJ3v0E4EiY0yxMaYFWALM6LiAMWalMebYVHhrgHTvxlSBpriyjqTYCOIjw6yOok7i2F9bu3TkTdDzpOjTgNIOj8vcz53MzcA7HR5HikieiKwRkSs6W0FE5ruXyausrPQgkvJ3JYd0jht/d+z7U6InZIOeV0/GishcIBd4uMPT/Y0xucB1wKMiMvDE9Ywxi4wxucaY3OTkZG9GUhYpPlTPQC16vxYfGUZSbISekA0BnhR9OZDR4XG6+7njiMjFwM+A6caYL6+rNsaUu/9bDHwIjDuDvCoAVDW0cKS+hawkPRHr79pH3ugefbDzpOjXAdkiMkBEwoHZwHGjZ0RkHPAU7SVf0eH5HiIS4f48CTgP2Oqt8Mo/7dITsQFjYLKOpQ8FXRa9McYJ3AGsALYBS40xhSLyoIgcG0XzMBALvHLCMMphQJ6IbARWAr81xmjRB7ljhwJ0aKX/y0qK5Uh9C1UNLVZHUT7k8GQhY8zbwNsnPPeLDp9ffJL1PgNGnUlAFXiKD9XjsAnpPXQyM3937IK2XZX15PQPtziN8hW9MlZ5XXFlHf16RRNm139e/k5vKxga9CdRed2OA7UM6a33iQ0E/XpGE+GwsfNArdVRlA9p0Suvqm92sudIA0P7xFsdRXnAYbcxuHcc27Xog5oWvfKqnQdrMQaGpeoefaAYlhrH9gM1VsdQPqRFr7zq2J7hsFTdow8UQ/vEc6iuhcpava1gsNKiV161fX8NsREO0hJ1xE2gGOr+60v36oOXFr3yqm0HahnSJ05vHxhAjp1P2b5fj9MHKy165TXGGLbvr2FoHz0+H0h6xoTTOz6CbbpHH7S06JXX7KtuoqbJyVA9Ph9whvaJZ5vu0QctLXrlNdv3t+8RDtM9+oAzNDWOoopaWl16t6lgpEWvvObYiJshWvQBZ1ifeFpdRmeyDFJa9Mprtu2vIaNnFHF6V6mAc2w4rI68CU5a9Mprth+o1StiA1RWcgxhdtHj9EFKi155RVOri+LKOj0+H6DC7DYGpegVssFKi155RVFFHW0GHXETwIb1idOx9EFKi155xTb3iBsdQx+4hqbGcaCmiaP1ehOSYKNFr7xi+4FaIsNs9O+ltw8MVF9eIaszWQYdLXrlFVv31TC4dxx2nfogYB2b86ZwX7XFSZS3adGrM9bY4iJ/71Fy+/e0Ooo6AylxkfTvFc3qXYetjqK8TItenbE1JYdpcbYxZUiy1VHUGZo8OJnPdh2m2emyOoryIi16dcY+2lFJZJiNiQN0jz7QTR6cTGOri7zdR62OorxIi16dsY93VnJ2Vi8iw+xWR1Fn6OysXoTbbXy0s9LqKMqLtOjVGdl7uIHiQ/VMHqyHbYJBTISDCQN68OGOCqujKC/Soldn5KOd7YUwZUiKxUmUt0wZnMLOg3Xsq2q0OoryEi16dUY+2llJv57RZPaKtjqK8pLJ7pPqH+vhm6DhUdGLyFQR2SEiRSKysJPX7xWRrSKySUTeF5H+HV6bJyJfuD/meTO8slaz08Vnuw4zeXAyIjp+Plhkp8SSmhCpx+mDSJdFLyJ24HFgGjAcmCMiw09YbAOQa4wZDSwDfu9etydwP3AWMBG4X0R6eC++slL+7qM0tLj0+HyQEREmD05m1ReH9EYkQcKTPfqJQJExptgY0wIsAWZ0XMAYs9IY0+B+uAZId39+GfCuMeaIMeYo8C4w1TvRldU+2llJmF04Z2Avq6MoL5s8OJnaZicFpVVWR1Fe4EnRpwGlHR6XuZ87mZuBd05zXRUgjDG8v72C3P49iYlwWB1Hedm5g5Kw24T3th20OoryAq+ejBWRuUAu8PAprjdfRPJEJK+yUo8LBoINpVUUVdTx7TF9rY6ifCAhKozJg5N5fX05Tj18E/A8KfpyIKPD43T3c8cRkYuBnwHTjTHNp7KuMWaRMSbXGJObnKzHewPBkrV7iQ63M32sFn2wmj0hg4raZj7YrmPqA50nRb8OyBaRASISDswGlndcQETGAU/RXvId/1WsAC4VkR7uk7CXup9TAay2qZW3Nu7n26P7EquHbYLWhUNTSImLYMm60q4XVn6ty6I3xjiBO2gv6G3AUmNMoYg8KCLT3Ys9DMQCr4hIgYgsd697BHiI9l8W64AH3c+pALZ84z4aW13MnpjR9cIqYDnsNq7NzeDDHRXsr9aLpwKZR7tjxpi3gbdPeO4XHT6/+GvWfRp4+nQDKv+zZG0pQ/vEMTYj0eooysdmTcjgsZVFLF1Xxl0XZ1sdR50mvTJWnZIt5dVsLq9m9oQMvUgqBGT0jGZSdhJL80pxtRmr46jTpEWvTsmSdXuJcNi4clx61wuroDB7Qj/Kqxr55AsdEReotOiVx5paXby5YR+Xj0olITrM6jiqm1wyvDe9YsJ5WU/KBiwteuWxFYUHqG12MjNH9+ZDSbjDxoyxaby/rYKj9S1Wx1GnQYteeWxZfhlpiVGcnaVTHoSaq3PSaHG18damfVZHUadBi1555EB1E58WHeKq8WnYbHoSNtSM6JvAsNR4Xs0vszqKOg1a9Mojr20oo83A1eP1sE2ounp8GhvLqvniYK3VUdQp0qJXXTLG8Gp+Gbn9e5CZFGN1HGWRK8al4bAJy9brXn2g0aJXXdpYVs2uynqu0ZOwIS0pNoIpQ5J5Y0O5jqkPMFr0qkvL8kuJcNi4fHSq1VGUxa4en87BmmYdUx9gtOjV12pscfHWxv1MHdmH+EgdOx/qLhyWQmJ0GK/k6eGbQKJFr77WaxvKqG5s5bqJ/ayOovxAhMPOtbkZ/LvwAOVVOtFZoNCiVyfV1mZYvKqEUWkJTBzQ0+o4yk/MOzcTgGc+221pDuU5LXp1Uh/urKC4sp5bJg3QCczUl9ISo5g2sg8vfb6Xuman1XGUB7To1UktXlVCn/hILh+lJ2HV8W6ZlEVts5NX8nT+m0CgRa86tXVfDZ8WHWbeuZmE2fWfiTre2IxEcvr34OlPS3SoZQDQn2DVqac/LSEqzK4nYdVJ3XL+AEqPNPLu1oNWR1Fd0KJXX1FR08Tygn3MzE3X6YjVSV06og/pPaL4+yfFVkdRXdCiV1/x9Ke7cba1cfP5A6yOovyY3SbcfP4A8vYcJX+P3gran2nRq+PUNrXywpo9TBuVSv9eOq+N+nqzJmSQGB3Gkx/pXr0/06JXx3lp7V5qm50s+MZAq6OoABAd7uDGczJ5d+tBiirqrI6jTkKLXn2p2eli8aoSzhvUi1HpCVbHUQFi3jn9iQyzsejjXVZHUSehRa++9GbBPg7WNHOr7s2rU9ArNoJrczN4fUM5B2uarI6jOqFFr4D26Q4WfVzM8NR4JmUnWR1HBZhbzs/C1WZ4elWJ1VFUJ7ToFQCvri+jqKKOBVMG6nQH6pT16xXNt8f05dnVe3SyMz+kRa+obmjlt+9sZ3y/RL6l0x2o03TfZUMwGB56a6vVUdQJPCp6EZkqIjtEpEhEFnby+jdEZL2IOEXkmhNec4lIgftjubeCK+/5w392cLShhYeuGKk3/lanLb1HNHdemM2/Cw/w4Y4Kq+OoDrosehGxA48D04DhwBwRGX7CYnuB7wAvdvIWjcaYse6P6WeYV3nZ5rJqnv98Dzeek8mIvjrSRp2ZWyYNICsphgeWF9LsdFkdR7l5skc/ESgyxhQbY1qAJcCMjgsYY3YbYzYBbT7IqHykrc3wX29uoVdMBPdcMtjqOCoIRDjs/HLGCHYfbmCRXkTlNzwp+jSg41ykZe7nPBUpInkiskZEruhsARGZ714mr7JS70XZXZ7+tISC0ip+evlQEqJ0ThvlHZOyk/nmqFT+srKI7QdqrI6j6J6Tsf2NMbnAdcCjIvKVQdrGmEXGmFxjTG5ycnI3RFKF+6r5/b93cMnw3lw57lR+byvVtV/OGEF8ZBh3vVRAU6sewrGaJ0VfDmR0eJzufs4jxphy93+LgQ+BcaeQT/lAY4uLu5YUkBgdxu+uHq3DKZXXJcVG8IeZo9lxsJbfvrPd6jghz5OiXwdki8gAEQkHZgMejZ4RkR4iEuH+PAk4D9CxVxb777e3UVRRxx+vHUPPmHCr46ggNWVICt85N5N/fLablToKx1JdFr0xxgncAawAtgFLjTGFIvKgiEwHEJEJIlIGzASeEpFC9+rDgDwR2QisBH5rjNGit9BnRYd4bs0evjdpAJOy9TCZ8q2F04YypHcc972yiXq9v6xlxBj/ug1Ybm6uycvLszpGUDLGcPUTn7G/uomV/28KkWF2qyOpELB+71Gu+utn/GjqEG6fMsjqOEFLRPLd50O/Qq+MDSEf7axk/d4q7rhwkJa86jbj+/XggiHJLPq4mNqmVqvjhCQt+hBhjOGRd3eSlhjFzJyMrldQyovuuWQwVQ2t/OPT3VZHCUla9CHig+0VbCyr5gcXDSLcod921b1Gpydy8bDe/O2TYqobda++u+lPfAgwxvDIezvp1zOaq8anWx1Hhai7L86mpsmpUxlbQIs+BLy+oZwt5TX84KJswuz6LVfWGJmWwNQRfVi8qoTSIw1Wxwkp+lMf5Lbuq+Gnr29mQmYPrhjb1+o4KsQtnDYUARY8n69XzHYjLfogVtXQwq3P55EQFcbj14/HoXvzymKZSTE8Onsshe4dEH8b3h2s9Cc/SLnaDD9YUsDB6maemJtDSlyk1ZGUAuCiYb255+LBvLa+nGc+2211nJCgRR+k/vCfHXy8s5JfzhjB+H49rI6j1HHuvHAQFw/rza/+tY3Piw9bHSfoadEHoX9t2s8TH+5izsR+zJnYz+o4Sn2FzSb8adYY+vWK5vYX1rNP7zPrU1r0QWbHgVruW7aR8f0SeWD6iTcCU8p/xEeGseiGXJqdbdymJ2d9Sos+iFQ3tDL/uTxiIxw8OTeHCIdOc6D826CUWB6ZNZaNZdX8/I0tenLWR7Tog4Qxhh++spF9VY08MXc8KfF68lUFhkuG9+aui7JZll/Gy+tKu15BnTIt+iDx7Oo9vLftIAunDSOnf0+r4yh1Su66KJvzByXxwFuFFFXUWh0n6GjRB4Gt+2r49dvbuHBoCjedl2l1HKVOmc0m/OnaMcSEO7jjxQ16vN7LtOgDXEOLkztfWk9iVBgPX6O3BVSBKyU+kj/MHMP2A7X85u1tVscJKlr0Ae6hf26j+FA9j8waS6/YCKvjKHVGLhiaws3nD+CZ1Xt4b+tBq+MEDS36APb+toO8tHYv8ydlcd6gJKvjKOUVP5o6hKF94lj42iYO1zVbHScoaNEHqMN1zfz41c0M7RPHvZcOtjqOUl4T4bDz6Oyx1DQ6dT4cL9GiD0DGGH72+hZqGlt5ZNZYHS+vgs7QPvH88NLBrCg8yKvry62OE/C06APQa+vL+XfhAX546WCGpcZbHUcpn7hlUhYTB/TkgeWFOn/9GdKiDzC7D9Xzize3MHFAT26ZlGV1HKV8xm4T/jhzDAB3v1yA09VmcaLApUUfQFqcbfxgyQYcdhuPzhqL3aZDKVVwy+gZza+vHEn+nqP8z/tfWB0nYGnRB5CHV2xnU1k1v79mNH0To6yOo1S3mDE2jZk56Ty2sojPdh2yOk5A0qIPEB/uqOBvn5Rww9n9uWxEH6vjKNWtfjljBAOSYrjn5QKO1LdYHSfgeFT0IjJVRHaISJGILOzk9W+IyHoRcYrINSe8Nk9EvnB/zPNW8FCyr6qRe5duZGifOH72zWFWx1Gq20WHO/jLnHEcrW/l7pcLcLXpkMtT0WXRi4gdeByYBgwH5ojIiROd7wW+A7x4wro9gfuBs4CJwP0iorc7OgXNThe3v7CeFmcbj18/nsgwHUqpQtOIvgk8MH0EH++s5M96vP6UeLJHPxEoMsYUG2NagCXAjI4LGGN2G2M2ASeeFr8MeNcYc8QYcxR4F5jqhdwh49f/2kZBaRUPXzOagcmxVsdRylJzJmZw9fh0/vzBF6zcUWF1nIDhSdGnAR0niS5zP+cJj9YVkfkikicieZWVlR6+dfB7Y0M5z67ew/cmDWDaqFSr4yhlORHhV1eMZEjvOO5eUqDj6z3kFydjjTGLjDG5xpjc5ORkq+P4hS3l1Sx8bRMTM3vyo6lDrY6jlN+ICrfz5Nwc2oxh/nP5NLQ4rY7k9zwp+nIgo8PjdPdznjiTdUNWZW0z33s2j57R4Tx+/XjC7H7x+1gpv5GZFMOf54xj+4Ea7ntlk86H0wVPGmQdkC0iA0QkHJgNLPfw/VcAl4pID/dJ2Evdz6mTaHa6uO35fI42tLDoxlyS43TqYaU6c8GQFBZOHcq/Nu/nsQ+KrI7j17osemOME7iD9oLeBiw1xhSKyIMiMh1ARCaISBkwE3hKRArd6x4BHqL9l8U64EH3c6oTbW2Gn7++hbw9R3n4mjGMTEuwOpJSfm3+N7K4clwaf3x3J29v3m91HL8l/vYnT25ursnLy7M6RrdztRkWvrqJV/LL+MFF2dx7iU49rJQnmlpdXPe3NWwsq+bRWWP59pi+VkeyhIjkG2NyO3tND/76gVZXG3e/XMAr+WXcdVE291ycbXUkpQJGZJidZ26aSE7/Hty1ZANL80q7XinEaNFbrKHFyW3Pr+etjftYOG0o91wyWO/7qtQpiosM45nvTuS8QUn8aNkmFq8q0RO0HWjRW2h/dSMzn1zN+9sP8uCMESyYPNDqSEoFrKhwO3+fl8tlI3rz0D+38rM3ttCqUxsDWvSWKSitYvpjn7LncAOL5+Vy4zmZVkdSKuBFOOw8cX0OCyYP5MXP9zLv6bVUNegkaFr0Fvhg+0FmPbWayDAbr91+LhcO7W11JKWChs0mLJw2lD/OHEPe7qNc9cRn7KtqtDqWpbTou9mbBeXMfzaf7N6xvH77eQzuHWd1JKWC0tU56Tx380Qqa5q55onPKK6sszqSZbTou9Fza/Zw98sF5PTvwUvfO5ukWL0YSilfOiurFy/NP5tmZxszn1zNlvJqqyNZQou+mzy+soj/emMLFw1N4ZmbJhIXGWZ1JKVCwsi0BF5ZcA4RDhtz/raGvN2hd82mFr2PGWP47TvbeXjFDq4Y25cn5ubonPJKdbOs5FiW3XYuybER3LB4LR/vDK1ZcrXofaitzfCLNwt58qNdXH9WP/507VidoEwpi/RNjOLlW88hMymGW57J499bQmfKBG0dH2lrM/z09c08t2YPCyYP5FdXjMRm0wuhlLJSclwES753NiPT4vn+ixtCZn4cLXofaGsz/OyNzSxZV8qdFw7ix1OH6NWuSvmJhOgwnr35LMZlJHLnSxt4JwTKXovey9pLfgsvrS3ljgsGca9OaaCU34mNcPCPmyYy1l32wX4YR4veyx55bycvrd3L7VMG8sNLteSV8lexEQ7+8d0JjE5P4M6XNrAuiEfjaNF70ZsF5fzlgyJm5WZw32V6uEYpfxcXGcb/fmciGT2iufW5/KC9B60WvZds2HuU+5ZtYuKAnjx0xUgteaUCREJ0GH+fl4vT1cYtz+RR29RqdSSv06L3gv3Vjcx/Lp8+8ZE8OTeHcIf+b1UqkGQlx/LE3ByKKuu4a0kBrrbgmuJYG+kMNTtdLHh+PY0tLhbPy6VnTLjVkZRSp+G8QUk88O3hfLC9gkff22l1HK/Soj9D979ZyMbSKv547RiydYIypQLa3LP7c21uOn/5oIgVhQesjuM1WvRn4MXP97JkXfswystG9LE6jlLqDIkID84YyZj0BH64dCNFFcEx46UW/WnK33OU+5dvYfLgZO7RG3krFTQiw+w8MTeHCIeNW5/LoyYITs5q0Z+GsqMN3PpcHn0To/if2WOx69QGSgWVvolRPHbdePYcbuDOFzfgDPBbEmrRn6L6Zie3PJNHs7ONxfMmkBitJ1+VCkbnDOzFgzNG8tHOSv777e1WxzkjDqsDBBJXm+GuJQV8UVHHP747gUEpsVZHUkr50HVn9aOooo6nPy1hUEos153Vz+pIp0WL3kPGGH79r228t+0gD84YwaTsZKsjKaW6wU8vH0rxoTp+8eYW0npEMXlw4P3se3ToRkSmisgOESkSkYWdvB4hIi+7X/9cRDLdz2eKSKOIFLg/nvRy/m7z1MfFPP1pCd89L5Mbz8m0Oo5Sqps47Db+Mmcc2b3juO35fApKq6yOdMq6LHoRsQOPA9OA4cAcERl+wmI3A0eNMYOAR4DfdXhtlzFmrPtjgZdyd6tX8kr57TvbmT6mL//1zRM3XSkV7OIiw3jmuxPoFRvOTf9Yx64Au9G4J3v0E4EiY0yxMaYFWALMOGGZGcAz7s+XARdJkEz28p/CAyx8bTPnD0riDzPH6M1DlApRKfGRPHvTWQhw4+K17KtqtDqSxzwp+jSgtMPjMvdznS5jjHEC1UAv92sDRGSDiHwkIpM6+wIiMl9E8kQkr7LSf+7l+GZBObe/sJ6RfeN58gadw0apUDcgKYZ/fHciNY2tzHxyNSWH6q2O5BFfN9d+oJ8xZhxwL/CiiMSfuJAxZpExJtcYk5uc7B8nOl78fC93v1xATv8ePH/LWcRG6HlrpRSMSk/gpfln09jqYuaTq9m2v8bqSF3ypOjLgYwOj9Pdz3W6jIg4gATgsDGm2RhzGMAYkw/sAvz6MlJjDI+vLOKnr2/mgiEpPHPTROIiw6yOpZTyIyPTElh66zk4bMKsp1bzefFhqyN9LU+Kfh2QLSIDRCQcmA0sP2GZ5cA89+fXAB8YY4yIJLtP5iIiWUA2UOyd6N7X1Ori7pcLeHjFDq4Y25cn5+YQGWa3OpZSyg8NSonllQXnkBQXwfV//5yX1u61OtJJdVn07mPudwArgG3AUmNMoYg8KCLT3YstBnqJSBHth2iODcH8BrBJRApoP0m7wBjjl/frOljTxKxFa3izYB/3XTaER2aN1WPySqmvldEzmtdvP49zByXxk9c288DyQr+cLkGM8a8J9nNzc01eXl63fs11u49w+wvrqW928qdrxzJ1pM5EqZTynNPVxm/e2c7iVSWcndWTx64bT1JsRLdmEJF8Y0xuZ6+F9C6rMYb//bSEOYvWEBNu57Xbz9WSV0qdMofdxn99azh/nDmGDXur+NafV7F+71GrY30pZIveGMNPXtvML9/aygVDU1h+5/kM7fOVAUFKKeWxq3PSee32cwlztJ+k/dem/VZHAkK46B9fWcSSdaXcPmUgT83NIV5H1iilvGBE3wT+ecckxqQncs/SAjb4wZ59SBb9O5v384f/7OTKcWncd9kQvdpVKeVVCdFhLLoxl97xEcx/Lt/yq2hDrui3lFdzz9ICxvdL5DdXjSJIZmpQSvmZnjHhPD1vAk0tLm55Jo+GFqdlWUKq6JtaXdz2Qj69YiJ46oZcHSOvlPKp7N5x/Pm6cWw/UMOv/rXNshwhVfTPrt5N6ZFGfn/NaJLjunfok1IqNF0wJIUbz8lkydq9fHGw1pIMIVP0VQ0tPPZBEVOGJHPeoCSr4yilQsgPLsomJtzBb9+x5paEIVP0f/mgiLpmJz+ZNszqKEqpENMzJpzbLxjE+9srWL2r++fFCYmi33u4gWdX72ZmTgZD+sRZHUcpFYK+e14mfRMi+e+3t9HW1r0zEoRE0T/8nx3YbcI9l/j1xJlKqSAWGWbnh5cOYXN5NW9t2tetXzvoi35/dSP/3LSP75w7gD4JkVbHUUqFsCvHpTG4dyyLV5V069cN+qJ/fUM5xsDsCRldL6yUUj5kswnX5mawqay6W0fgBHXRG2N4Nb+M3P49yEyKsTqOUkoxY2wadpuwbH1Zt33NoC76jWXV7Kqs55qcdKujKKUUAMlxEVwwJJk3NpTj6qaTskFd9K/mlxHhsHH56FSroyil1JeuHp/OwZpmVhUd6pavF7RF3+x0sXzjPi4b0UdnplRK+ZULh6WQEBXGq/ndc/gmaIv+g20VVDe2crUetlFK+ZkIh53pY/qyovAANU2tPv96QVv0y/LL6B0fwfk63YFSyg9dnZNOs7OtW25OEpRFX3a0gQ93VnLFuPaz20op5W/GpCcwKCWWFz7fg6/v3R2URb/o42JsAt85N9PqKEop1SkR4ZbzB7ClvIZPvvDtSdmgK/qK2iaWrCvl6vHppCZEWR1HKaVO6qrx6aQmRPL4yiKffp2gK/rFq0pwutq4dfJAq6MopdTXCnfY+N6kLD4vOULe7iM++zpBVfTVDa08v3oP3xzdlwF6JaxSKgDMnphBz5hwn+7VB1XRP7N6N/UtLm6fonvzSqnAEB3u4ObzB7ByRyWF+6p98jU8KnoRmSoiO0SkSEQWdvJ6hIi87H79cxHJ7PDaT9zP7xCRy7yY/Tj1zU6e/rSEi4amMCw13ldfRimlvG7u2f2Ji3Dw15W7fPL+jq4WEBE78DhwCVAGrBOR5caYrR0Wuxk4aowZJCKzgd8Bs0RkODAbGAH0Bd4TkcHGGJe3N6Su2cm5A3tx8/lZ3n5rpZTyqYSoMBZMGUhTqwtjDCLeHRbeZdEDE4EiY0wxgIgsAWYAHYt+BvCA+/NlwGPSnnQGsMQY0wyUiEiR+/1Weyf+/+kdH8lfr8/x9tsqpVS3+P4Fg3z23p4cukkDSjs8LnM/1+kyxhgnUA308nBdpZRSPuQXJ2NFZL6I5IlIXmVlpdVxlFIqqHhS9OVAx9szpbuf63QZEXEACcBhD9fFGLPIGJNrjMlNTk72PL1SSqkueVL064BsERkgIuG0n1xdfsIyy4F57s+vAT4w7ZM3LAdmu0flDACygbXeia6UUsoTXZ6MNcY4ReQOYAVgB542xhSKyINAnjFmObAYeM59svUI7b8McC+3lPYTt07g+74YcaOUUurkxNezpp2q3Nxck5eXZ3UMpZQKKCKSb4zJ7ew1vzgZq5RSyne06JVSKsj53aEbEakE9pzBWyQB3XPHXd8Klu0A3RZ/FSzbEizbAWe2Lf2NMZ0OW/S7oj9TIpJ3suNUgSRYtgN0W/xVsGxLsGwH+G5b9NCNUkoFOS16pZQKcsFY9IusDuAlwbIdoNvir4JlW4JlO8BH2xJ0x+iVUkodLxj36JVSSnWgRa+UUkEu6IpeRB4SkU0iUiAi/xGRvlZnOl0i8rCIbHdvz+sikmh1ptMlIjNFpFBE2kQk4IbCdXU7zUAiIk+LSIWIbLE6y5kQkQwRWSkiW93/tu6yOtPpEpFIEVkrIhvd2/JLr75/sB2jF5F4Y0yN+/MfAMONMQssjnVaRORS2mcCdYrI7wCMMT+2ONZpEZFhQBvwFPD/jDEBM6GR+3aaO+lwO01gzgm30wwYIvINoA541hgz0uo8p0tEUoFUY8x6EYkD8oErAvH74r4jX4wxpk5EwoBVwF3GmDXeeP+g26M/VvJuMUDA/iYzxvzHfccugDW0z+cfkIwx24wxO6zOcZq+vJ2mMaYFOHY7zYBkjPmY9llmA5oxZr8xZr3781pgGwF6BzvTrs79MMz94bXuCrqiBxCRX4tIKXA98Aur83jJTcA7VocIUXpLTD8nIpnAOOBzi6OcNhGxi0gBUAG8a4zx2rYEZNGLyHsisqWTjxkAxpifGWMygBeAO6xN+/W62hb3Mj+jfT7/F6xL2jVPtkUpbxORWOBV4O4T/qIPKMYYlzFmLO1/uU8UEa8dVuvyxiP+yBhzsYeLvgC8DdzvwzhnpKttEZHvAN8CLjJ+fkLlFL4vgcajW2Kq7uc+nv0q8IIx5jWr83iDMaZKRFYCUwGvnDAPyD36ryMi2R0ezgC2W5XlTInIVOBHwHRjTIPVeUKYJ7fTVN3MfQJzMbDNGPMnq/OcCRFJPjaqTkSiaD/x77XuCsZRN68CQ2gf4bEHWGCMCci9L/etGSNov9E6wJoAHkF0JfAXIBmoAgqMMZdZGuoUiMjlwKP83+00f21totMnIi8BU2ifEvcgcL8xZrGloU6DiJwPfAJspv3nHeCnxpi3rUt1ekRkNPAM7f++bMBSY8yDXnv/YCt6pZRSxwu6QzdKKaWOp0WvlFJBToteKaWCnBa9UkoFOS16pZQKclr0SikV5LTolVIqyP1/XD+okqjobWoAAAAASUVORK5CYII=\n",
107+ "text/plain": [
108+ "<Figure size 432x288 with 1 Axes>"
109+ ]
110+ },
111+ "metadata": {
112+ "needs_background": "light"
113+ },
114+ "output_type": "display_data"
115+ }
116+ ],
117+ "source": [
118+ "plt.plot(e0,dos)"
119+ ]
120+ },
121+ {
122+ "cell_type": "code",
123+ "execution_count": 8,
124+ "id": "ec75ef83",
125+ "metadata": {},
126+ "outputs": [],
127+ "source": [
128+ "ef, wght, iteration=libtetrabz.fermieng(bvec, eig, 0.4)"
129+ ]
130+ },
131+ {
132+ "cell_type": "code",
133+ "execution_count": 9,
134+ "id": "e88213bb",
135+ "metadata": {},
136+ "outputs": [
137+ {
138+ "data": {
139+ "text/plain": [
140+ "0.3999999999670125"
141+ ]
142+ },
143+ "execution_count": 9,
144+ "metadata": {},
145+ "output_type": "execute_result"
146+ }
147+ ],
148+ "source": [
149+ "wght.sum()"
150+ ]
151+ },
152+ {
153+ "cell_type": "code",
154+ "execution_count": 10,
155+ "id": "4c259f40",
156+ "metadata": {},
157+ "outputs": [
158+ {
159+ "data": {
160+ "text/plain": [
161+ "-0.3498372919857502"
162+ ]
163+ },
164+ "execution_count": 10,
165+ "metadata": {},
166+ "output_type": "execute_result"
167+ }
168+ ],
169+ "source": [
170+ "ef"
171+ ]
172+ },
173+ {
174+ "cell_type": "code",
175+ "execution_count": 11,
176+ "id": "09446911",
177+ "metadata": {},
178+ "outputs": [
179+ {
180+ "data": {
181+ "text/plain": [
182+ "28"
183+ ]
184+ },
185+ "execution_count": 11,
186+ "metadata": {},
187+ "output_type": "execute_result"
188+ }
189+ ],
190+ "source": [
191+ "iteration"
192+ ]
193+ },
194+ {
195+ "cell_type": "code",
196+ "execution_count": null,
197+ "id": "e28c3ef4",
198+ "metadata": {},
199+ "outputs": [],
200+ "source": [
201+ "with open(\"sc.frmsf\", \"w\") as f:\n",
202+ " print(ng[0], ng[1], ng[2], file=f)\n",
203+ " print(1, file=f)\n",
204+ " print(nb, file=f)\n",
205+ " for ii in range(3):\n",
206+ " print(bvec[ii, 0], bvec[ii, 1], bvec[ii, 2], file=f)\n",
207+ " for ib in range(nb):\n",
208+ " for i0 in range(ng[0]):\n",
209+ " for i1 in range(ng[1]):\n",
210+ " for i2 in range(ng[2]):\n",
211+ " print(eig[i0, i1, i2, ib]-ef, file=f)"
212+ ]
213+ },
214+ {
215+ "cell_type": "code",
216+ "execution_count": null,
217+ "id": "ba33dd26",
218+ "metadata": {},
219+ "outputs": [],
220+ "source": [
221+ "# Lindhard function"
222+ ]
223+ },
224+ {
225+ "cell_type": "code",
226+ "execution_count": null,
227+ "id": "11206e36",
228+ "metadata": {},
229+ "outputs": [],
230+ "source": []
231+ }
232+ ],
233+ "metadata": {
234+ "kernelspec": {
235+ "display_name": "Python 3 (ipykernel)",
236+ "language": "python",
237+ "name": "python3"
238+ },
239+ "language_info": {
240+ "codemirror_mode": {
241+ "name": "ipython",
242+ "version": 3
243+ },
244+ "file_extension": ".py",
245+ "mimetype": "text/x-python",
246+ "name": "python",
247+ "nbconvert_exporter": "python",
248+ "pygments_lexer": "ipython3",
249+ "version": "3.7.8"
250+ }
251+ },
252+ "nbformat": 4,
253+ "nbformat_minor": 5
254+}
Show on old repository browser