• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revisión8cae309368cc0f56524723509d3a9b1455193f97 (tree)
Tiempo2021-11-02 01:34:28
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Bug fix

Cambiar Resumen

Diferencia incremental

--- a/README.md
+++ b/README.md
@@ -1,5 +1,16 @@
11 # libtetrabz
22
3-This is a simple example package. You can use
4-[Github-flavored Markdown](https://guides.github.com/features/mastering-markdown/)
5-to write your content.
3+Libtetrabz is a library which parform efficiently the Brillouin-zone
4+integration in the electronic structure calculation in a solid by
5+using the tetrahedron method.
6+
7+## For citing libtetrabz
8+
9+We would appreciate if you cite the following article in your
10+research with libtetrabz.
11+
12+"Improved tetrahedron method for the Brillouin-zone integration applicable to response functions",
13+
14+M. Kawamura, Y. Gohda, and S. Tsuneyuki, Phys. Rev. B 89, 094515 (2014).
15+
16+https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.094515
--- a/src/libtetrabz/__init__.py
+++ b/src/libtetrabz/__init__.py
@@ -20,12 +20,12 @@
2020 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
23-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
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/libtetrabz/libtetrabz_common.py
+++ b/src/libtetrabz/libtetrabz_common.py
@@ -132,15 +132,15 @@ def libtetrabz_initialize(ng, bvec):
132132 ikv = numpy.empty([nk*6, 20, 3], dtype=numpy.int_)
133133 ikv0 = numpy.empty(3, dtype=numpy.int_)
134134 nt = 0
135- for i3 in range(ng[2]):
136- for i2 in range(ng[1]):
137- for i1 in range(ng[0]):
135+ for i2 in range(ng[2]):
136+ for i1 in range(ng[1]):
137+ for i0 in range(ng[0]):
138138 #
139139 for it in range(6):
140140 #
141141 for ii in range(20):
142142 #
143- ikv0[0:3] = numpy.array([i1, i2, i3]) + ivvec[it, ii, 0:3]
143+ ikv0[0:3] = [i0, i1, i2] + ivvec[it, ii, 0:3]
144144 ikv[nt, ii, 0:3] = ikv0[0:3] % ng[0:3]
145145 #
146146 nt += 1
@@ -150,16 +150,16 @@ def libtetrabz_initialize(ng, bvec):
150150 def libtetrabz_tsmall_a1(e):
151151 """Cut small tetrahedron A1
152152 """
153- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
154- for ii in range(4):
155- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
153+ a10 = (0.0 - e[0]) / (e[1] - e[0])
154+ a20 = (0.0 - e[0]) / (e[2] - e[0])
155+ a30 = (0.0 - e[0]) / (e[3] - e[0])
156156 #
157- v = a[1, 0] * a[2, 0] * a[3, 0]
157+ v = a10 * a20 * a30
158158 #
159- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
160- [a[0, 1], a[1, 0], 0.0, 0.0],
161- [a[0, 2], 0.0, a[2, 0], 0.0],
162- [a[0, 3], 0.0, 0.0, a[3, 0]]])
159+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
160+ [1.0 - a10, a10, 0.0, 0.0],
161+ [1.0 - a20, 0.0, a20, 0.0],
162+ [1.0 - a30, 0.0, 0.0, a30]])
163163 return v, tsmall
164164
165165
@@ -171,15 +171,15 @@ def libtetrabz_tsmall_b1(e):
171171 """
172172 #
173173 #
174- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
175- for ii in range(4):
176- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
174+ a13 = (0.0 - e[3]) / (e[1] - e[3])
175+ a20 = (0.0 - e[0]) / (e[2] - e[0])
176+ a30 = (0.0 - e[0]) / (e[3] - e[0])
177177 #
178- v = a[2, 0] * a[3, 0] * a[1, 3]
179- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
180- [a[0, 2], 0.0, a[2, 0], 0.0],
181- [a[0, 3], 0.0, 0.0, a[3, 0]],
182- [0.0, a[1, 3], 0.0, a[3, 1]]])
178+ v = a20 * a30 * a13
179+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
180+ [1.0 - a20, 0.0, a20, 0.0],
181+ [1.0 - a30, 0.0, 0.0, a30],
182+ [0.0, a13, 0.0, 1.0 - a13]])
183183 return v, tsmall
184184
185185
@@ -189,16 +189,15 @@ def libtetrabz_tsmall_b2(e):
189189 :param e:
190190 :return:
191191 """
192- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
193- for ii in range(4):
194- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
192+ a21 = (0.0 - e[1]) / (e[2] - e[1])
193+ a31 = (0.0 - e[1]) / (e[3] - e[1])
195194 #
196- v = a[2, 1] * a[3, 1]
195+ v = a21 * a31
197196 #
198- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
199- [0.0, 1.0, 0.0, 0.0],
200- [0.0, a[1, 2], a[2, 1], 0.0],
201- [0.0, a[1, 3], 0.0, a[3, 1]]])
197+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
198+ [0.0, 1.0, 0.0, 0.0],
199+ [0.0, 1.0 - a21, a21, 0.0],
200+ [0.0, 1.0 - a31, 0.0, a31]])
202201 return v, tsmall
203202
204203
@@ -208,16 +207,16 @@ def libtetrabz_tsmall_b3(e):
208207 :param e:
209208 :return:
210209 """
211- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
212- for ii in range(4):
213- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
210+ a12 = (0.0 - e[2]) / (e[1] - e[2])
211+ a20 = (0.0 - e[0]) / (e[2] - e[0])
212+ a31 = (0.0 - e[1]) / (e[3] - e[1])
214213 #
215- v = a[1, 2] * a[2, 0] * a[3, 1]
214+ v = a12 * a20 * a31
216215 #
217- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
218- [a[0, 2], 0.0, a[2, 0], 0.0],
219- [0.0, a[1, 2], a[2, 1], 0.0],
220- [0.0, a[1, 3], 0.0, a[3, 1]]])
216+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
217+ [1.0-a20, 0.0, a20, 0.0],
218+ [0.0, a12, 1.0-a12, 0.0],
219+ [0.0, 1.0-a31, 0.0, a31]])
221220 return v, tsmall
222221
223222
@@ -227,16 +226,14 @@ def libtetrabz_tsmall_c1(e):
227226 :param e:
228227 :return:
229228 """
230- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
231- for ii in range(4):
232- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
229+ a32 = (0.0 - e[2]) / (e[3] - e[2])
233230 #
234- v = a[3, 2]
231+ v = a32
235232 #
236- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
237- [0.0, 1.0, 0.0, 0.0],
238- [0.0, 0.0, 1.0, 0.0],
239- [0.0, 0.0, a[2, 3], a[3, 2]]])
233+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
234+ [0.0, 1.0, 0.0, 0.0],
235+ [0.0, 0.0, 1.0, 0.0],
236+ [0.0, 0.0, 1.0 - a32, a32]])
240237 return v, tsmall
241238
242239
@@ -246,16 +243,15 @@ def libtetrabz_tsmall_c2(e):
246243 :param e:
247244 :return:
248245 """
249- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
250- for ii in range(4):
251- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
246+ a23 = (0.0 - e[3]) / (e[2] - e[3])
247+ a31 = (0.0 - e[1]) / (e[3] - e[1])
252248 #
253- v = a[2, 3] * a[3, 1]
249+ v = a23 * a31
254250 #
255- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
256- [0.0, 1.0, 0.0, 0.0],
257- [0.0, a[1, 3], 0.0, a[3, 1]],
258- [0.0, 0.0, a[2, 3], a[3, 2]]])
251+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
252+ [0.0, 1.0, 0.0, 0.0],
253+ [0.0, 1.0 - a31, 0.0, a31],
254+ [0.0, 0.0, a23, 1.0 - a23]])
259255 return v, tsmall
260256
261257
@@ -265,16 +261,16 @@ def libtetrabz_tsmall_c3(e):
265261 :param e:
266262 :return:
267263 """
268- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
269- for ii in range(4):
270- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
264+ a23 = (0.0 - e[3]) / (e[2] - e[3])
265+ a13 = (0.0 - e[3]) / (e[1] - e[3])
266+ a30 = (0.0 - e[0]) / (e[3] - e[0])
271267 #
272- v = a[2, 3] * a[1, 3] * a[3, 0]
268+ v = a23 * a13 * a30
273269 #
274- tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
275- [a[0, 3], 0.0, 0.0, a[3, 0]],
276- [0.0, a[1, 3], 0.0, a[3, 1]],
277- [0.0, 0.0, a[2, 3], a[3, 2]]])
270+ tsmall = numpy.array([[1.0, 0.0, 0.0, 0.0],
271+ [1.0 - a30, 0.0, 0.0, a30],
272+ [0.0, a13, 0.0, 1.0 - a13],
273+ [0.0, 0.0, a23, 1.0 - a23]])
278274 return v, tsmall
279275
280276
@@ -284,16 +280,16 @@ def libtetrabz_triangle_a1(e):
284280 :param e:
285281 :return:
286282 """
287- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
288- for ii in range(4):
289- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
283+ a10 = (0.0 - e[0]) / (e[1] - e[0])
284+ a20 = (0.0 - e[0]) / (e[2] - e[0])
285+ a30 = (0.0 - e[0]) / (e[3] - e[0])
290286 #
291- # v = 3.0 * a[1,0] * a[2,0] * a[3,0] / (0.0 - e(1))
292- v = 3.0 * a[1, 0] * a[2, 0] / (e(4) - e(1))
287+ # v = 3.0 * a[1,0] * a[2,0] * a[3,0] / (0.0 - e[0])
288+ v = 3.0 * a10 * a20 / (e[3] - e[0])
293289 #
294- tsmall = numpy.array([[a[0, 1], a[1, 0], 0.0, 0.0],
295- [a[0, 2], 0.0, a[2, 0], 0.0],
296- [a[0, 3], 0.0, 0.0, a[3, 0]]])
290+ tsmall = numpy.array([[1.0 - a10, a10, 0.0, 0.0],
291+ [1.0 - a20, 0.0, a20, 0.0],
292+ [1.0 - a30, 0.0, 0.0, a30]])
297293 return v, tsmall
298294
299295
@@ -303,16 +299,16 @@ def libtetrabz_triangle_b1(e):
303299 :param e:
304300 :return:
305301 """
306- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
307- for ii in range(4):
308- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
302+ a30 = (0.0 - e[0]) / (e[3] - e[0])
303+ a13 = (0.0 - e[3]) / (e[1] - e[3])
304+ a20 = (0.0 - e[0]) / (e[2] - e[0])
309305 #
310- # v = 3.0 * a[2,0] * a[3,0] * a[1,3] / (0.0 - e(1))
311- v = 3.0 * a[3, 0] * a[1, 3] / (e(3) - e(1))
306+ # v = 3.0 * a[2,0] * a[3,0] * a[1,3] / (0.0 - e[0])
307+ v = 3.0 * a30 * a13 / (e[2] - e[0])
312308 #
313- tsmall = numpy.array([[a[0, 2], 0.0, a[2, 0], 0.0],
314- [a[0, 3], 0.0, 0.0, a[3, 0]],
315- [0.0, a[1, 3], 0.0, a[3, 1]]])
309+ tsmall = numpy.array([[1.0 - a20, 0.0, a20, 0.0],
310+ [1.0 - a30, 0.0, 0.0, a30],
311+ [0.0, a13, 0.0, 1.0 - a13]])
316312 return v, tsmall
317313
318314
@@ -321,16 +317,16 @@ def libtetrabz_triangle_b2(e):
321317 :param e:
322318 :return:
323319 """
324- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
325- for ii in range(4):
326- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
320+ a12 = (0.0 - e[2]) / (e[1] - e[2])
321+ a31 = (0.0 - e[1]) / (e[3] - e[1])
322+ a20 = (0.0 - e[0]) / (e[2] - e[0])
327323 #
328- # v = 3.0 * a[1,2] * a[2,0] * a[3,1] / (0.0 - e(1))
329- v = 3.0 * a[1, 2] * a[3, 1] / (e[2] - e[0])
324+ # v = 3.0 * a[1,2] * a[2,0] * a[3,1] / (0.0 - e[0])
325+ v = 3.0 * a12 * a31 / (e[2] - e[0])
330326 #
331- tsmall = numpy.array([[a[0, 2], 0.0, a[2, 0], 0.0],
332- [0.0, a[1, 2], a[2, 1], 0.0],
333- [0.0, a[1, 3], 0.0, a[3, 1]]])
327+ tsmall = numpy.array([[1.0 - a20, 0.0, a20, 0.0],
328+ [0.0, a12, 1.0 - a12, 0.0],
329+ [0.0, 1.0 - a31, 0.0, a31]])
334330 return v, tsmall
335331
336332
@@ -339,14 +335,14 @@ def libtetrabz_triangle_c1(e):
339335 :param e:
340336 :return:
341337 """
342- a = numpy.empty(shape=(4, 4), dtype=numpy.float_)
343- for ii in range(4):
344- a[0:4, ii] = (0.0 - e[ii]) / (e[0:4] - e[ii])
338+ a03 = (0.0 - e[3]) / (e[0] - e[3])
339+ a13 = (0.0 - e[3]) / (e[1] - e[3])
340+ a23 = (0.0 - e[3]) / (e[2] - e[3])
345341 #
346- # v = 3.0 * a[0,3] * a[1,3] * a[2,3] / (e(4) - 0.0)
347- v = 3.0 * a[0, 3] * a[1, 3] / (e[3] - e[2])
342+ # v = 3.0 * a[0,3] * a[1,3] * a[2,3] / (e[3] - 0.0)
343+ v = 3.0 * a03 * a13 / (e[3] - e[2])
348344 #
349- tsmall = numpy.array([[a[0, 3], 0.0, 0.0, a[3, 0]],
350- [0.0, a[1, 3], 0.0, a[3, 1]],
351- [0.0, 0.0, a[2, 3], a[3, 2]]])
345+ tsmall = numpy.array([[a03, 0.0, 0.0, 1.0 - a03],
346+ [0.0, a13, 0.0, 1.0 - a13],
347+ [0.0, 0.0, a23, 1.0 - a23]])
352348 return v, tsmall
--- a/src/libtetrabz/libtetrabz_dbldelta.py
+++ b/src/libtetrabz/libtetrabz_dbldelta.py
@@ -21,7 +21,7 @@
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
2323 import numpy
24-import libtetrabz_common
24+from . import libtetrabz_common
2525
2626
2727 def dbldelta(bvec, eig1, eig2):
@@ -37,23 +37,22 @@ def dbldelta(bvec, eig1, eig2):
3737 nk = ng.prod(0)
3838 nb = eig1.shape[3]
3939 wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
40-
40+ #
4141 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
42-
43- eig1t = numpy.empty([20, nb], dtype=numpy.float_)
44- eig2t = numpy.empty([20, nb], dtype=numpy.float_)
45-
42+ #
4643 for it in range(6 * nk):
4744 #
45+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
46+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
4847 for ii in range(20):
4948 eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
5049 eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
51-
5250 ei1 = wlsm.dot(eig1t)
5351 ej1 = wlsm.dot(eig2t)
52+ #
53+ w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
5454 for ib in range(nb):
5555 #
56- w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
5756 e = ei1[0:4, ib]
5857 indx = e.argsort(0)
5958 e.sort(0)
@@ -66,7 +65,7 @@ def dbldelta(bvec, eig1, eig2):
6665 #
6766 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
6867 w2 = libtetrabz_dbldelta2(ej2)
69- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
68+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
7069 #
7170 elif e[1] < 0.0 <= e[2]:
7271 #
@@ -75,14 +74,14 @@ def dbldelta(bvec, eig1, eig2):
7574 if v > thr:
7675 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7776 w2 = libtetrabz_dbldelta2(ej2)
78- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
77+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
7978 #
8079 v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e)
8180 #
8281 if v > thr:
8382 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8483 w2 = libtetrabz_dbldelta2(ej2)
85- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
84+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
8685 #
8786 elif e[2] < 0.0 < e[3]:
8887 #
@@ -91,15 +90,16 @@ def dbldelta(bvec, eig1, eig2):
9190 if v > thr:
9291 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9392 w2 = libtetrabz_dbldelta2(ej2)
94- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
93+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
9594 #
9695 else:
9796 continue
98- #
99- for ii in range(20):
100- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
97+ #
98+ for ii in range(20):
99+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb] += w1.dot(wlsm[:, ii])
101100 #
102101 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
102+ return wght
103103
104104
105105 def libtetrabz_dbldelta2(ej):
--- a/src/libtetrabz/libtetrabz_dblstep.py
+++ b/src/libtetrabz/libtetrabz_dblstep.py
@@ -21,7 +21,7 @@
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
2323 import numpy
24-import libtetrabz_common
24+from . import libtetrabz_common
2525
2626
2727 def dblstep(bvec, eig1, eig2):
@@ -39,23 +39,21 @@ def dblstep(bvec, eig1, eig2):
3939
4040 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
4141
42- eig1t = numpy.empty([20, nb], dtype=numpy.float_)
43- eig2t = numpy.empty([20, nb], dtype=numpy.float_)
44-
4542 thr = 1.0e-8
4643 #
4744 for it in range(6 * nk):
4845 #
46+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
47+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
4948 for ii in range(20):
5049 eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
5150 eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
52-
5351 ei1 = wlsm.dot(eig1t)
5452 ej1 = wlsm.dot(eig2t)
5553 #
54+ w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
5655 for ib in range(nb):
5756 #
58- w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
5957 e = ei1[0:4, ib]
6058 indx = e.argsort(0)
6159 e.sort(0)
@@ -65,84 +63,76 @@ def dblstep(bvec, eig1, eig2):
6563 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
6664 #
6765 if v > thr:
68- #
6966 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7067 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7168 w2 = libtetrabz_dblstep2(ei2, ej2)
72- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
69+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
7370 #
7471 elif e(2) <= 0.0 < e(3):
7572 #
7673 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
7774 #
7875 if v > thr:
79- #
8076 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8177 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8278 w2 = libtetrabz_dblstep2(ei2, ej2)
83- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
79+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
8480 #
8581 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8682 #
8783 if v > thr:
88- #
8984 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9085 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9186 w2 = libtetrabz_dblstep2(ei2, ej2)
92- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
87+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
9388 #
9489 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9590 #
9691 if v > thr:
97- #
9892 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9993 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
10094 w2 = libtetrabz_dblstep2(ei2, ej2)
101- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
95+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
10296 #
10397 elif e(3) <= 0.0 < e(4):
10498 #
10599 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
106100 #
107101 if v > thr:
108- #
109102 ei2 = tsmall.dot(ei1[indx[0:4], ib])
110103 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
111104 w2 = libtetrabz_dblstep2(ei2, ej2)
112- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
105+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
113106 #
114107 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
115108 #
116109 if v > thr:
117- #
118110 ei2 = tsmall.dot(ei1[indx[0:4], ib])
119111 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
120112 w2 = libtetrabz_dblstep2(ei2, ej2)
121- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
113+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
122114 #
123115 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
124116 #
125117 if v > thr:
126- #
127118 ei2 = tsmall.dot(ei1[indx[0:4], ib])
128119 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
129120 w2 = libtetrabz_dblstep2(ei2, ej2)
130- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
121+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
131122 #
132123 elif e(4) <= 0.0:
133- #
134124 ei2 = ei1[0:4, ib]
135125 ej2 = ej1[0:4, 0:nb]
136126 w2 = libtetrabz_dblstep2(ei2, ej2)
137- w1[0:nb, 0:4] += w2[0:nb, 0:4]
138- #
127+ w1[ib, 0:nb, 0:4] += w2[0:nb, 0:4]
139128 else:
140129 continue
141130 #
142- for ii in range(20):
143- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
144-
131+ for ii in range(20):
132+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb] += w1.dot(wlsm[:, ii])
133+ #
145134 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
135+ return wght
146136
147137
148138 def libtetrabz_dblstep2(ei1, ej1):
--- a/src/libtetrabz/libtetrabz_dos.py
+++ b/src/libtetrabz/libtetrabz_dos.py
@@ -21,16 +21,16 @@
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
2323 import numpy
24-import libtetrabz_common
24+from . import libtetrabz_common
2525
2626
2727 def dos(bvec, eig, e0):
2828 """
2929 Compute Dos : Delta(E - E1)
30- :param bvec:
31- :param eig:
32- :param e0:
33- :return:
30+ :param ndarray([3, 3], float) bvec: Reciprocal-lattice vectors
31+ :param ndarray([ng0, ng1, ng2, nb], float) eig: Energy eigenvalue
32+ :param ndarray(ne, float) e0: Energy grid
33+ :return ndarray([ng0, ng1, ng2, nb, ne], float) wght: Weight
3434 """
3535
3636 ng = numpy.array(eig.shape[0:3])
@@ -41,47 +41,49 @@ def dos(bvec, eig, e0):
4141
4242 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, ne], dtype=numpy.float_)
4343
44- eigt = numpy.empty([20, nb], dtype=numpy.float_)
4544 for it in range(6*nk):
4645 #
46+ eigt = numpy.empty([20, nb], dtype=numpy.float_)
4747 for ii in range(20):
4848 eigt[ii, 0:nb] = eig[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
49-
5049 ei1 = wlsm.dot(eigt)
5150 #
51+ w1 = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
52+ #
5253 for ib in range(nb):
5354 #
54- w1 = numpy.zeros([ne, 4], dtype=numpy.float_)
5555 e = ei1[0:4, ib]
5656 indx = e.argsort(0)
57- e.sort(0)
57+ e = e[indx[0:4]]
5858 #
5959 for ie in range(ne):
6060 #
6161 if e[0] < e0[ie] <= e[1] or e[0] <= e0[ie] < e[1]:
6262 #
6363 v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e[0:4] - e0[ie])
64- w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
64+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
6565 #
6666 elif e[1] < e0[ie] <= e[2] or e[1] <= e0[ie] < e[2]:
6767 #
6868 v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e[0:4] - e0[ie])
69- w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
69+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
7070 #
7171 v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e[0:4] - e0[ie])
72- w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
72+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
7373 #
7474 elif e[2] < e0[ie] <= e[3] or e[2] <= e0[ie] < e[3]:
7575 #
7676 v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e[0:4] - e0[ie])
77- w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
77+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
7878 #
7979 else:
8080 continue
8181 #
82- for ii in range(20):
83- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:ne] += w1.dot(wlsm)
84- #
82+ #
83+ #
84+ for ii in range(20):
85+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:ne] += w1.dot(wlsm[:, ii])
86+ #
8587 #
8688 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:ne] /= (6.0 * nk)
8789 return wght
@@ -90,31 +92,30 @@ def dos(bvec, eig, e0):
9092 def intdos(bvec, eig, e0):
9193 """
9294 Compute integrated Dos : theta(E - E1)
93- :param bvec:
94- :param eig:
95- :param e0:
96- :return:
95+ :param ndarray([3, 3], float) bvec: Reciprocal-lattice vectors
96+ :param ndarray([ng0, ng1, ng2, nb], float) eig: Energy eigenvalue
97+ :param ndarray(ne, float) e0: Energy grid
98+ :return ndarray([ng0, ng1, ng2, nb, ne], float) wght: Weight
9799 """
98100 ng = numpy.array(eig.shape[0:3])
99101 nk = ng.prod(0)
100102 nb = eig.shape[3]
101103 ne = e0.shape[0]
102104 wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
103-
104105 #
105106 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, ne], dtype=numpy.float_)
106107 #
107- eigt = numpy.empty([20, nb], dtype=numpy.float_)
108108 for it in range(6*nk):
109109 #
110+ eigt = numpy.empty([20, nb], dtype=numpy.float_)
110111 for ii in range(20):
111112 eigt[ii, 0:nb] = eig[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
112-
113113 ei1 = wlsm.dot(eigt)
114114 #
115+ w1 = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
116+ #
115117 for ib in range(nb):
116118 #
117- w1 = numpy.zeros([ne, 4], dtype=numpy.float_)
118119 e = ei1[0:4, ib]
119120 indx = e.argsort(0)
120121 e.sort(0)
@@ -124,40 +125,39 @@ def intdos(bvec, eig, e0):
124125 if e[0] <= e0[ie] < e[1] or e[0] < e0[ie] <= e[1]:
125126 #
126127 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e - e0[ie])
127- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
128+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
128129 #
129130 elif e[1] <= e0[ie] < e[2] or e[1] < e0[ie] <= e[2]:
130131 #
131132 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e - e0[ie])
132- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
133+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
133134 #
134135 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e - e0[ie])
135- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
136+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
136137 #
137138 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e - e0[ie])
138- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
139- #
139+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
140+ #
140141 elif e[2] <= e0[ie] < e[3] or e[2] < e0[ie] .AND. e0[ie] <= e[3]:
141142 #
142143 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e - e0[ie])
143- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
144+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
144145 #
145146 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e - e0[ie])
146- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
147+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
147148 #
148149 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e - e0[ie])
149- w1[ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
150- #
151- elif e[3] <= e0(ie):
152- #
153- w1[ie, 0:4] = 0.25
150+ w1[ib, ie, indx[0:4]] += v * tsmall.sum(0) * 0.25
154151 #
152+ elif e[3] <= e0(ie):
153+ w1[ib, ie, 0:4] = 0.25
155154 else:
156- #
157155 continue
158156 #
159- for ii in range(20):
160- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:ne] += w1.dot(wlsm)
161-
157+ #
158+ #
159+ for ii in range(20):
160+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:ne] += w1.dot(wlsm[:, ii])
161+ #
162162 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:ne] /= (6.0 * nk)
163163 return wght
--- a/src/libtetrabz/libtetrabz_fermigr.py
+++ b/src/libtetrabz/libtetrabz_fermigr.py
@@ -21,7 +21,7 @@
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
2323 import numpy
24-import libtetrabz_common
24+from . import libtetrabz_common
2525
2626
2727 def fermigr(bvec, eig1, eig2, e0):
@@ -41,22 +41,21 @@ def fermigr(bvec, eig1, eig2, e0):
4141
4242 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.float_)
4343
44- eig1t = numpy.empty([20, nb], dtype=numpy.float_)
45- eig2t = numpy.empty([20, nb], dtype=numpy.float_)
4644 thr = 1.0e-10
4745 #
4846 for it in range(6 * nk):
4947 #
48+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
49+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
5050 for ii in range(20):
5151 eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
5252 eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
53-
5453 ei1 = wlsm.dot(eig1t)
5554 ej1 = wlsm.dot(eig2t)
5655 #
56+ w1 = numpy.zeros([nb, nb, ne, 4], dtype=numpy.float_)
5757 for ib in range(nb):
5858 #
59- w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
6059 e = ei1[0:4, ib]
6160 indx = e.argsort(0)
6261 e.sort(0)
@@ -70,7 +69,7 @@ def fermigr(bvec, eig1, eig2, e0):
7069 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7170 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7271 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
73- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
72+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
7473 #
7574 elif e[1] <= 0.0 < e[2]:
7675 #
@@ -81,7 +80,7 @@ def fermigr(bvec, eig1, eig2, e0):
8180 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8281 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8382 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
84- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
83+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
8584 #
8685 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8786 #
@@ -90,7 +89,7 @@ def fermigr(bvec, eig1, eig2, e0):
9089 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9190 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9291 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
93- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
92+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
9493 #
9594 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9695 #
@@ -99,7 +98,7 @@ def fermigr(bvec, eig1, eig2, e0):
9998 ei2 = tsmall.dot(ei1[indx[0:4], ib])
10099 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
101100 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
102- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
101+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
103102 #
104103 elif e[2] <= 0.0 < e[3]:
105104 #
@@ -110,7 +109,7 @@ def fermigr(bvec, eig1, eig2, e0):
110109 ei2 = tsmall.dot(ei1[indx[0:4], ib])
111110 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
112111 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
113- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
112+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
114113 #
115114 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
116115 #
@@ -119,7 +118,7 @@ def fermigr(bvec, eig1, eig2, e0):
119118 ei2 = tsmall.dot(ei1[indx[0:4], ib])
120119 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
121120 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
122- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
121+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
123122 #
124123 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
125124 #
@@ -128,22 +127,23 @@ def fermigr(bvec, eig1, eig2, e0):
128127 ei2 = tsmall.dot(ei1[indx[0:4], ib])
129128 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
130129 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
131- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
130+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
132131 #
133132 elif e[3] <= 0.0:
134133 #
135134 ei2 = ei1[0:4, ib]
136135 ej2 = ej1[0:4, 0:nb]
137136 w2 = libtetrabz_fermigr2(e0, ei2, ej2)
138- w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
137+ w1[ib, 0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
139138 #
140139 else:
141140 continue
142141 #
143- for ii in range(20):
144- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
145- #
142+ for ii in range(20):
143+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb, 0:ne] += w1.dot(wlsm[:, ii])
144+ #
146145 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
146+ return wght
147147
148148
149149 def libtetrabz_fermigr2(e0, ei1, ej1):
--- a/src/libtetrabz/libtetrabz_occ.py
+++ b/src/libtetrabz/libtetrabz_occ.py
@@ -21,15 +21,15 @@
2121 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 #
2323 import numpy
24-import libtetrabz_common
24+from . import libtetrabz_common
2525
2626
2727 def fermieng(bvec, eig, nelec):
2828 """
2929 Calculate Fermi energy
30- :param bvec:
31- :param eig:
32- :param nelec:
30+ :param ndarray([3, 3], float) bvec: Reciprocal-lattice vectors
31+ :param ndarray([ng0, ng1, ng2, nb], float) eig: Energy eigenvalue
32+ :param float nelec:
3333 :return:
3434 """
3535 #
@@ -52,8 +52,8 @@ def fermieng(bvec, eig, nelec):
5252 #
5353 # convergence check
5454 #
55- if sumkmid - nelec < eps:
56- return ef, wght
55+ if abs(sumkmid - nelec) < eps:
56+ return ef, wght, iteration
5757 elif sumkmid < nelec:
5858 elw = ef
5959 else:
@@ -62,72 +62,70 @@ def fermieng(bvec, eig, nelec):
6262 raise ValueError("libtetrabz_fermieng")
6363
6464
65-def occ(bvec, eig):
65+def occ(bvec=numpy.array([1.0, 0.0, 0.0]), eig=numpy.array([0.0])):
6666 """
6767 Main SUBROUTINE for occupation : Theta(EF - E1)
68- :param bvec:
69- :param eig:
70- :return:
68+ :param ndarray([3, 3], float) bvec: Reciprocal-lattice vectors
69+ :param ndarray([ng0, ng1, ng2, nb], float) eig: Energy eigenvalue
70+ :return ndarray([ng0, ng1, ng2, nb], float) wght: Weight
7171 """
7272 #
7373 ng = numpy.array(eig.shape[0:3])
7474 nk = ng.prod(0)
7575 nb = eig.shape[3]
7676 wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
77-
77+ #
7878 wght = numpy.zeros([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)
79-
80- eigt = numpy.empty([20, nb], dtype=numpy.float_)
79+ #
8180 for it in range(6 * nk):
8281 #
82+ eigt = numpy.empty([20, nb], dtype=numpy.float_)
8383 for ii in range(20):
8484 eigt[ii, 0:nb] = eig[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
85-
8685 ei1 = wlsm.dot(eigt)
8786 #
87+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
8888 for ib in range(nb):
8989 #
90- w1 = numpy.zeros(4, dtype=numpy.float_)
9190 e = ei1[0:4, ib]
9291 indx = e.argsort(0)
9392 e.sort(0)
9493 #
95- if e(1) <= 0.0 < e(2):
94+ if e[0] <= 0.0 < e[1]:
9695 #
9796 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
98- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
97+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
9998 #
100- elif e(2) <= 0.0 < e(3):
99+ elif e[1] <= 0.0 < e[2]:
101100 #
102101 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
103- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
102+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
104103 #
105104 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
106- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
105+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
107106 #
108107 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
109- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
108+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
110109 #
111- elif e(3) <= 0.0 < e(4):
110+ elif e[2] <= 0.0 < e[3]:
112111 #
113112 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
114- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
113+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
115114 #
116115 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
117- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
116+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
118117 #
119118 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
120- w1[indx[0:4]] += v * tsmall.sum(0) * 0.25
121- #
122- elif e(4) <= 0.0:
123- #
124- w1[0:4] = 0.25
119+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
125120 #
121+ elif e[3] <= 0.0:
122+ w1[ib, 0:4] = 0.25
126123 else:
127124 continue
128125 #
129- for ii in range(20):
130- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib] += w1.dot(wlsm)
131- #
126+ #
127+ for ii in range(20):
128+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb] += w1.dot(wlsm[:, ii])
129+ #
132130 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb] /= (6.0 * nk)
133131 return wght
--- a/src/libtetrabz/libtetrabz_polcmplx.py
+++ b/src/libtetrabz/libtetrabz_polcmplx.py
@@ -22,7 +22,7 @@
2222 #
2323 import math
2424 import numpy
25-import libtetrabz_common
25+from . import libtetrabz_common
2626
2727
2828 def polcmplx(bvec, eig1, eig2, e0):
@@ -41,24 +41,21 @@ def polcmplx(bvec, eig1, eig2, e0):
4141 wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
4242
4343 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.complex_)
44-
45- eig1t = numpy.empty([20, nb], dtype=numpy.float_)
46- eig2t = numpy.empty([20, nb], dtype=numpy.float_)
47-
4844 #
4945 thr = 1.0e-8
5046 for it in range(6 * nk):
5147 #
48+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
49+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
5250 for ii in range(20):
5351 eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
5452 eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
55-
5653 ei1 = wlsm.dot(eig1t)
5754 ej1 = wlsm.dot(eig2t)
5855 #
56+ w1 = numpy.zeros([nb, nb, ne, 4], dtype=numpy.float_)
5957 for ib in range(nb):
6058 #
61- w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
6259 e = ei1[0:4, ib]
6360 indx = e.argsort(0)
6461 e.sort(0)
@@ -68,91 +65,78 @@ def polcmplx(bvec, eig1, eig2, e0):
6865 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
6966 #
7067 if v > thr:
71- #
7268 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7369 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7470 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
75- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
76- #
71+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
7772 #
7873 elif e(2) <= 0.0 < e[2]:
7974 #
8075 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
8176 #
8277 if v > thr:
83- #
8478 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8579 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8680 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
87- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
88- #
81+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
8982 #
9083 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
9184 #
9285 if v > thr:
93- #
9486 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9587 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9688 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
97- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
98- #
89+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
9990 #
10091 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
10192 #
10293 if v > thr:
103- #
10494 ei2 = tsmall.dot(ei1[indx[0:4], ib])
10595 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
10696 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
107- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
108- #
97+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
10998 #
11099 elif e[2] <= 0.0 < e[3]:
111100 #
112101 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
113102 #
114103 if v > thr:
115- #
116104 ei2 = tsmall.dot(ei1[indx[0:4], ib])
117105 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
118106 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
119- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
120- #
107+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
121108 #
122109 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
123110 #
124111 if v > thr:
125- #
126112 ei2 = tsmall.dot(ei1[indx[0:4], ib])
127113 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
128114 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
129- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
130- #
115+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
131116 #
132117 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
133118 #
134119 if v > thr:
135- #
136120 ei2 = tsmall.dot(ei1[indx[0:4], ib])
137121 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
138122 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
139- w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
140- #
123+ w1[ib, 0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
141124 #
142125 elif e[3] <= 0.0:
143126 #
144127 ei2 = ei1[0:4, ib]
145128 ej2 = ej1[0:4, 0:nb]
146129 w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
147- w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
130+ w1[ib, 0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
148131 #
149132 else:
150133 continue
151134 #
152- for ii in range(20):
153- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
154- #
135+ for ii in range(20):
136+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb, 0:ne] += w1.dot(wlsm[:, ii])
137+ #
155138 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
139+ return wght
156140
157141
158142 def libtetrabz_polcmplx2(e0, ei1, ej1):
--- a/src/libtetrabz/libtetrabz_polstat.py
+++ b/src/libtetrabz/libtetrabz_polstat.py
@@ -22,7 +22,7 @@
2222 #
2323 import math
2424 import numpy
25-import libtetrabz_common
25+from . import libtetrabz_common
2626
2727
2828 def polstat(bvec, eig1, eig2):
@@ -42,106 +42,97 @@ def polstat(bvec, eig1, eig2):
4242
4343 wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
4444
45- eig1t = numpy.empty([20, nb], dtype=numpy.float_)
46- eig2t = numpy.empty([20, nb], dtype=numpy.float_)
4745 for it in range(6 * nk):
4846 #
47+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
48+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
4949 for ii in range(20):
5050 eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
5151 eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
52-
5352 ei1 = wlsm.dot(eig1t)
5453 ej1 = wlsm.dot(eig2t)
5554 #
55+ w1 = numpy.zeros([nb, nb, 4], dtype=numpy.float_)
5656 for ib in range(nb):
5757 #
58- w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
5958 e = ei1[0:4, ib]
6059 indx = e.argsort(0)
6160 e.sort(0)
62-
61+ #
6362 if e[0] <= 0.0 < e[1]:
6463 #
6564 v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
6665 #
6766 if v > thr:
68- #
6967 ei2 = tsmall.dot(ei1[indx[0:4], ib])
7068 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
7169 w2 = libtetrabz_polstat2(ei2, ej2)
72- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
70+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
7371 #
7472 elif e[1] <= 0.0 < e[2]:
7573 #
7674 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
7775 #
7876 if v > thr:
79- #
8077 ei2 = tsmall.dot(ei1[indx[0:4], ib])
8178 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
8279 w2 = libtetrabz_polstat2(ei2, ej2)
83- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
80+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
8481 #
8582 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
8683 #
8784 if v > thr:
88- #
8985 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9086 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
9187 w2 = libtetrabz_polstat2(ei2, ej2)
92- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
88+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
9389 #
9490 v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
9591 #
9692 if v > thr:
97- #
9893 ei2 = tsmall.dot(ei1[indx[0:4], ib])
9994 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
10095 w2 = libtetrabz_polstat2(ei2, ej2)
101- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
96+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
10297 #
10398 elif e[2] <= 0.0 < e[3]:
10499 #
105100 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
106101 #
107102 if v > thr:
108- #
109103 ei2 = tsmall.dot(ei1[indx[0:4], ib])
110104 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
111105 w2 = libtetrabz_polstat2(ei2, ej2)
112- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
106+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
113107 #
114108 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
115109 #
116110 if v > thr:
117- #
118111 ei2 = tsmall.dot(ei1[indx[0:4], ib])
119112 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
120113 w2 = libtetrabz_polstat2(ei2, ej2)
121- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
114+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
122115 #
123116 v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
124117 #
125118 if v > thr:
126- #
127119 ei2 = tsmall.dot(ei1[indx[0:4], ib])
128120 ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
129121 w2 = libtetrabz_polstat2(ei2, ej2)
130- w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
131- #
132- elif e[3] <= 0.0:
122+ w1[ib, 0:nb, indx[0:4]] += v * w2.dot(tsmall)
133123 #
124+ elif e[3] <= 0.0:
134125 ei2 = ei1[0:4, ib]
135126 ej2 = ej1[0:4, 0:nb]
136127 w2 = libtetrabz_polstat2(ei2, ej2)
137- w1[0:nb, 0:4] += w2[0:nb, 0:4]
138- #
128+ w1[ib, 0:nb, 0:4] += w2[0:nb, 0:4]
139129 else:
140130 continue
141131 #
142- for ii in range(20):
143- wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
144- #
132+ #
133+ for ii in range(20):
134+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb, 0:nb] += w1.dot(wlsm[:, ii])
135+ #
145136 wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
146137
147138
Show on old repository browser