OSDN Git Service

Port remaining functions
authorMitsuaki Kawamura <kawamitsuaki@gmail.com>
Thu, 21 Oct 2021 16:55:50 +0000 (01:55 +0900)
committerMitsuaki Kawamura <kawamitsuaki@gmail.com>
Thu, 21 Oct 2021 16:55:50 +0000 (01:55 +0900)
LICENSE
README.md
pyproject.toml
src/libtetrabz/__init__.py
src/libtetrabz/libtetrabz_dbldelta.py [new file with mode: 0644]
src/libtetrabz/libtetrabz_dblstep.py [new file with mode: 0644]
src/libtetrabz/libtetrabz_dos.py [moved from src/libtetrabz/libtetrabz_dos_mod.py with 98% similarity]
src/libtetrabz/libtetrabz_fermigr.py [new file with mode: 0644]
src/libtetrabz/libtetrabz_occ.py [moved from src/libtetrabz/libtetrabz_occ_mod.py with 97% similarity]
src/libtetrabz/libtetrabz_polcmplx.py [new file with mode: 0644]
src/libtetrabz/libtetrabz_polstat.py [moved from src/libtetrabz/libtetrabz_polstat_mod.py with 99% similarity]

diff --git a/LICENSE b/LICENSE
index 96f1555..4c48664 100644 (file)
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2018 The Python Packaging Authority
+Copyright (c) Copyright (c) 2014 Mitsuaki Kawamura
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
index ad9ff9a..021ceae 100644 (file)
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Example Package
+# libtetrabz
 
 This is a simple example package. You can use
 [Github-flavored Markdown](https://guides.github.com/features/mastering-markdown/)
index 374b58c..a3163eb 100644 (file)
@@ -2,5 +2,6 @@
 requires = [
     "setuptools>=42",
     "wheel"
+    "numpy"
 ]
 build-backend = "setuptools.build_meta"
index e69de29..67cbdf7 100644 (file)
@@ -0,0 +1,31 @@
+#\r
+# Copyright (c) 2014 Mitsuaki Kawamura\r
+#\r
+# Permission is hereby granted, free of charge, to any person obtaining a\r
+# copy of this software and associated documentation files (the\r
+# "Software"), to deal in the Software without restriction, including\r
+# without limitation the rights to use, copy, modify, merge, publish,\r
+# distribute, sublicense, and/or sell copies of the Software, and to\r
+# permit persons to whom the Software is furnished to do so, subject to\r
+# the following conditions:\r
+#\r
+# The above copyright notice and this permission notice shall be included\r
+# in all copies or substantial portions of the Software.\r
+#\r
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS\r
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.\r
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY\r
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,\r
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE\r
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
+#\r
+from libtetrabz_dos import dos\r
+from libtetrabz_dos import intdos\r
+from libtetrabz_occ import occ\r
+from libtetrabz_occ import fermieng\r
+from libtetrabz_polstat import polstat\r
+from libtetrabz_fermigr import fermigr\r
+from libtetrabz_dbldelta import dbldelta\r
+from libtetrabz_dblstep import dblstep\r
+from libtetrabz_polcmplx import polcmplx\r
diff --git a/src/libtetrabz/libtetrabz_dbldelta.py b/src/libtetrabz/libtetrabz_dbldelta.py
new file mode 100644 (file)
index 0000000..d3b1435
--- /dev/null
@@ -0,0 +1,144 @@
+#
+# Copyright (c) 2014 Mitsuaki Kawamura
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+import numpy
+import libtetrabz_common
+
+
+def dbldelta(bvec, eig1, eig2):
+    """
+    Main SUBROUTINE for Delta(E1) * Delta(E2)
+    :param bvec:
+    :param eig1:
+    :param eig2:
+    :return:
+    """
+    thr = 1.0e-10
+    ng = numpy.array(eig1.shape[0:3])
+    nk = ng.prod(0)
+    nb = eig1.shape[3]
+    wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
+
+    wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
+
+    eig1t = numpy.empty([20, nb], dtype=numpy.float_)
+    eig2t = numpy.empty([20, nb], dtype=numpy.float_)
+
+    for it in range(6 * nk):
+        #
+        for ii in range(20):
+            eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+            eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+
+        ei1 = wlsm.dot(eig1t)
+        ej1 = wlsm.dot(eig2t)
+        for ib in range(nb):
+            #
+            w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
+            e = ei1[0:4, ib]
+            indx = e.argsort(0)
+            e.sort(0)
+            #
+            if e[0] < 0.0 <= e[1]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e)
+                #
+                if v > thr:
+                    #
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dbldelta2(ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e[1] < 0.0 <= e[2]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e)
+                #
+                if v > thr:
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dbldelta2(ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e)
+                #
+                if v > thr:
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dbldelta2(ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e[2] < 0.0 < e[3]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e)
+                #
+                if v > thr:
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dbldelta2(ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            else:
+                continue
+            #
+            for ii in range(20):
+                wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
+    #
+    wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
+
+
+def libtetrabz_dbldelta2(ej):
+    """
+    2nd step of tetrahedra method.
+    :param ej:
+    :return:
+    """
+    nb = ej.shape[1]
+    w = numpy.zeros([nb, 3], dtype=numpy.float_)
+    a = numpy.empty([3, 3], dtype=numpy.float_)
+    for ib in range(nb):
+        #
+        if abs(ej[0:3, ib]).max() < 1.0e-10:
+            raise ValueError("Nesting ##")
+        #
+        e = ej[0:3, ib]
+        indx = e.argsort(0)
+        e.sort(0)
+        #
+        for ii in range(3):
+            a[0:3, ii] = (0.0 - e[ii]) / (e[0:3] - e[ii])
+        #
+        if e[0] < 0.0 <= e[1] or e[0] <= 0.0 < e[1]:
+            #
+            # v = a[1, 0] * a[2, 0] / (0.0 - e[0])
+            v = a[1, 0] / (e[2] - e[0])
+            #
+            w[ib, indx[0]] = v * (a[0, 1] + a[0, 2])
+            w[ib, indx[1]] = v * a[1, 0]
+            w[ib, indx[2]] = v * a[2, 0]
+            #
+        elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
+            #
+            # v = a[1, 2] * a[1, 2] / (e[2] - 0.0)
+            v = a[1, 2] / (e[2] - e[0])
+            #
+            w[ib, indx[0]] = v * a[0, 2]
+            w[ib, indx[1]] = v * a[1, 2]
+            w[ib, indx[2]] = v * (a[2, 0] + a[2, 1])
+        #
+    return w
diff --git a/src/libtetrabz/libtetrabz_dblstep.py b/src/libtetrabz/libtetrabz_dblstep.py
new file mode 100644 (file)
index 0000000..2cd53cf
--- /dev/null
@@ -0,0 +1,203 @@
+#
+# Copyright (c) 2014 Mitsuaki Kawamura
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+import numpy
+import libtetrabz_common
+
+
+def dblstep(bvec, eig1, eig2):
+    """
+    Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
+    :param bvec:
+    :param eig1:
+    :param eig2:
+    :return:
+    """
+    ng = numpy.array(eig1.shape[0:3])
+    nk = ng.prod(0)
+    nb = eig1.shape[3]
+    wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
+
+    wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
+
+    eig1t = numpy.empty([20, nb], dtype=numpy.float_)
+    eig2t = numpy.empty([20, nb], dtype=numpy.float_)
+
+    thr = 1.0e-8
+    #
+    for it in range(6 * nk):
+        #
+        for ii in range(20):
+            eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+            eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+
+        ei1 = wlsm.dot(eig1t)
+        ej1 = wlsm.dot(eig2t)
+        #
+        for ib in range(nb):
+            #
+            w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
+            e = ei1[0:4, ib]
+            indx = e.argsort(0)
+            e.sort(0)
+            #
+            if e(1) <= 0.0 < e(2):
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e(2) <= 0.0 < e(3):
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e(3) <= 0.0 < e(4):
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_dblstep2(ei2, ej2)
+                    w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e(4) <= 0.0:
+                #
+                ei2 = ei1[0:4, ib]
+                ej2 = ej1[0:4, 0:nb]
+                w2 = libtetrabz_dblstep2(ei2, ej2)
+                w1[0:nb, 0:4] += w2[0:nb, 0:4]
+                #
+            else:
+                continue
+            #
+            for ii in range(20):
+                wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
+
+    wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
+
+
+def libtetrabz_dblstep2(ei1, ej1):
+    """
+    Tetrahedra method for theta( - de)
+    :param ei1:
+    :param ej1:
+    :return:
+    """
+    thr = 1.0e-8
+    nb = ej1.shape[1]
+    w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
+
+    for ib in range(nb):
+        #
+        e = - ei1[0:4] + ej1[0:4, ib]
+        indx = e.argsort(0)
+        e.sort(0)
+        #
+        if abs(e(1)) < thr and abs(e(4)) < thr:
+            #
+            # Theta(0) = 0.5
+            #
+            v = 0.5
+            w1[ib, 0:4] += v * 0.25
+            #
+        elif e(1) <= 0.0 < e(2) or e(1) < 0.0 <= e(2):
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+        elif e(2) <= 0.0 < e(3) or e(2) < 0.0 <= e(3):
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+        elif e(3) <= 0.0 < e(4) or e(3) < 0.0 <= e(4):
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+            w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
+            #
+        elif e(4) <= 0.0:
+            #
+            w1[ib, 0:4] += 0.25
+            #
+    return w1
similarity index 98%
rename from src/libtetrabz/libtetrabz_dos_mod.py
rename to src/libtetrabz/libtetrabz_dos.py
index 80a980f..5d23554 100644 (file)
@@ -24,7 +24,7 @@ import numpy
 import libtetrabz_common
 
 
-def libtetrabz_dos(bvec, eig, e0):
+def dos(bvec, eig, e0):
     """
     Compute Dos : Delta(E - E1)
     :param bvec:
@@ -87,7 +87,7 @@ def libtetrabz_dos(bvec, eig, e0):
     return wght
 
 
-def libtetrabz_intdos(bvec, eig, e0):
+def intdos(bvec, eig, e0):
     """
     Compute integrated Dos : theta(E - E1)
     :param bvec:
diff --git a/src/libtetrabz/libtetrabz_fermigr.py b/src/libtetrabz/libtetrabz_fermigr.py
new file mode 100644 (file)
index 0000000..0fe4453
--- /dev/null
@@ -0,0 +1,275 @@
+#
+# Copyright (c) 2014 Mitsuaki Kawamura
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+import numpy
+import libtetrabz_common
+
+
+def fermigr(bvec, eig1, eig2, e0):
+    """
+    Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
+    :param bvec:
+    :param eig1:
+    :param eig2:
+    :param e0:
+    :return:
+    """
+    ng = numpy.array(eig1.shape[0:3])
+    nk = ng.prod(0)
+    nb = eig1.shape[3]
+    ne = e0.shape[0]
+    wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
+
+    wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.float_)
+
+    eig1t = numpy.empty([20, nb], dtype=numpy.float_)
+    eig2t = numpy.empty([20, nb], dtype=numpy.float_)
+    thr = 1.0e-10
+    #
+    for it in range(6 * nk):
+        #
+        for ii in range(20):
+            eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+            eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+
+        ei1 = wlsm.dot(eig1t)
+        ej1 = wlsm.dot(eig2t)
+        #
+        for ib in range(nb):
+            #
+            w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
+            e = ei1[0:4, ib]
+            indx = e.argsort(0)
+            e.sort(0)
+            #
+            if e[0] <= 0.0 < e[1]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e[1] <= 0.0 < e[2]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e[2] <= 0.0 < e[3]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            elif e[3] <= 0.0:
+                #
+                ei2 = ei1[0:4, ib]
+                ej2 = ej1[0:4, 0:nb]
+                w2 = libtetrabz_fermigr2(e0, ei2, ej2)
+                w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
+                #
+            else:
+                continue
+            #
+            for ii in range(20):
+                wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
+    #
+    wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
+
+
+def libtetrabz_fermigr2(e0, ei1, ej1):
+    """
+    Tetrahedra method for theta( - E2)
+    :param e0:
+    :param ei1:
+    :param ej1:
+    :return:
+    """
+    #
+    ne = e0.shape[0]
+    nb = ej1.shape[1]
+    w1 = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
+    thr = 1.0e-8
+    #
+    for ib in range(nb):
+        #
+        e = -ej1[0:4, ib]
+        indx = e.argsort(0)
+        e.sort(0)
+        #
+        if e[0] <= 0.0 < e[1] or e[0] < 0.0 <= e[1]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_fermigr3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[3] <= 0.0:
+            #
+            de = ej1[0:4, ib] - ei1[0:4]
+            w2 = libtetrabz_fermigr3(e0, de)
+            w1[ib, 0:ne, 0:4] += w2
+            #
+    return w1
+
+
+def libtetrabz_fermigr3(e0, de):
+    #
+    ne = e0.shape[0]
+    e = de[0:4]
+    indx = e.argsort(0)
+    e.sort(0)
+    w1 = numpy.zeros([ne, 4], dtype=numpy.float_)
+    #
+    for ie in range(ne):
+        #
+        if e[0] < e0[ie] .AND. e0[ie] <= e[1]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e[0:4] - e0[ie])
+            w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
+            #
+        elif e[1] < e0[ie] .AND. e0[ie] <= e[2]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e[0:4] - e0[ie])
+            w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e[0:4] - e0[ie])
+            w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
+            #
+        elif e[2] < e0[ie] .AND. e0[ie] < e[3]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e[0:4] - e0[ie])
+            w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
+            #
+    return w1
similarity index 97%
rename from src/libtetrabz/libtetrabz_occ_mod.py
rename to src/libtetrabz/libtetrabz_occ.py
index 650bcd5..2c145ca 100644 (file)
@@ -24,7 +24,7 @@ import numpy
 import libtetrabz_common
 
 
-def libtetrabz_fermieng(bvec, eig, nelec):
+def fermieng(bvec, eig, nelec):
     """
     Calculate Fermi energy
     :param bvec:
@@ -47,7 +47,7 @@ def libtetrabz_fermieng(bvec, eig, nelec):
         #
         # Calc. # of electrons
         #
-        wght = libtetrabz_occ(bvec, eig-ef)
+        wght = occ(bvec, eig-ef)
         sumkmid = wght.sum()
         #
         # convergence check
@@ -62,7 +62,7 @@ def libtetrabz_fermieng(bvec, eig, nelec):
     raise ValueError("libtetrabz_fermieng")
 
 
-def libtetrabz_occ(bvec, eig):
+def occ(bvec, eig):
     """
     Main SUBROUTINE for occupation : Theta(EF - E1)
     :param bvec:
diff --git a/src/libtetrabz/libtetrabz_polcmplx.py b/src/libtetrabz/libtetrabz_polcmplx.py
new file mode 100644 (file)
index 0000000..c3b4a56
--- /dev/null
@@ -0,0 +1,622 @@
+#
+# Copyright (c) 2014 Mitsuaki Kawamura
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+import math
+import numpy
+import libtetrabz_common
+
+
+def polcmplx(bvec, eig1, eig2, e0):
+    """
+    Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
+    :param bvec:
+    :param eig1:
+    :param eig2:
+    :param e0:
+    :return:
+    """
+    ng = numpy.array(eig1.shape[0:3])
+    nk = ng.prod(0)
+    nb = eig1.shape[3]
+    ne = e0.shape[0]
+    wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
+
+    wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.complex_)
+
+    eig1t = numpy.empty([20, nb], dtype=numpy.float_)
+    eig2t = numpy.empty([20, nb], dtype=numpy.float_)
+
+    #
+    thr = 1.0e-8
+    for it in range(6 * nk):
+        #
+        for ii in range(20):
+            eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+            eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
+
+        ei1 = wlsm.dot(eig1t)
+        ej1 = wlsm.dot(eig2t)
+        #
+        for ib in range(nb):
+            #
+            w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
+            e = ei1[0:4, ib]
+            indx = e.argsort(0)
+            e.sort(0)
+            #
+            if e[0] <= 0.0 < e(2):
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+            elif e(2) <= 0.0 < e[2]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+            elif e[2] <= 0.0 < e[3]:
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+                v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+                #
+                if v > thr:
+                    #
+                    ei2 = tsmall.dot(ei1[indx[0:4], ib])
+                    ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
+                    w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                    w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                    #
+                #
+            elif e[3] <= 0.0:
+                #
+                ei2 = ei1[0:4, ib]
+                ej2 = ej1[0:4, 0:nb]
+                w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
+                w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
+                #
+            else:
+                continue
+            #
+            for ii in range(20):
+                wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
+            #
+    wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
+
+
+def libtetrabz_polcmplx2(e0, ei1, ej1):
+    """
+    Tetrahedra method for theta( - E2)
+    :param e0:
+    :param ei1:
+    :param ej1:
+    :return:
+    """
+    ne = e0.shape[0]
+    nb = ej1.shape[1]
+    thr = 1.0e-8
+    w1 = numpy.zeros([nb, ne, 4], dtype=numpy.complex_)
+
+    for ib in range(nb):
+        #
+        e = -ej1[0:4, ib]
+        indx = e.argsort(0)
+        e.sort(0)
+        #
+        if e[0] <= 0.0 < e(2) or e[0] < 0.0 <= e[1]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+            v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
+            #
+            if v > thr:
+                #
+                de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
+                w2 = libtetrabz_polcmplx3(e0, de)
+                w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
+                #
+            #
+        elif e[3] <= 0.0:
+            #
+            de = ej1[0:4, ib] - ei1[0:4]
+            w2 = libtetrabz_polcmplx3(e0, de)
+            w1[ib, 0:ne, 0:4] += w2
+
+    return w1
+
+
+def libtetrabz_polcmplx3(e0, de):
+    """
+    Tetrahedron method for delta(om - ep + e)
+    :param e0:
+    :param de:
+    :return:
+    """
+    ne = e0.shape[0]
+    e = de[0:4]
+    indx = e.argsort(0)
+    e.sort(0)
+    w1 = numpy.zeros([ne, 4], dtype=numpy.complex_)
+    w2 = numpy.empty([2, 4], dtype=numpy.float_)
+    #
+    for ie in range(ne):
+        #
+        # I don't know which one is better.
+        # The former is more stable.
+        # The latter is more accurate ?
+        #
+        w1[ie, 0:4] = 0.25 / (de[0:4] + e0[ie])
+        #
+        continue
+        #
+        x = (e[0:4] + e0[ie].real) / e0[ie].imag
+        # thr = maxval(de(1:4)) * 1d-3
+        thr = max(1.0e-3,  x.max() * 1.0e-2)
+        #
+        if abs(x[3] - x(3)) < thr:
+            if abs(x[3] - x(2)) < thr:
+                if abs(x[3] - x(1)) < thr:
+                    #
+                    # e[3] = e[2] = e[1] = e[0]
+                    #
+                    w2[0, 3] = 0.25 * x[3] / (1.0 + x[3]**2)
+                    w2[1, 3] = 0.25 / (1.0 + x[3]**2)
+                    w2[0:2, 2] = w2[0:2, 3]
+                    w2[0:2, 1] = w2[0:2, 3]
+                    w2[0:2, 0] = w2[0:2, 3]
+                    #
+                else:
+                    #
+                    # e[3] = e[2] = e[1]
+                    #
+                    w2[0:2, 3] = libtetrabz_polcmplx_1211(x[3], x(1))
+                    w2[0:2, 2] = w2[0:2, 3]
+                    w2[0:2, 1] = w2[0:2, 3]
+                    w2[0:2, 0] = libtetrabz_polcmplx_1222(x[0], x(4))
+                    #
+                    # IF(ANY(w2(1:2,1:4) < 0.0)):
+                    #   WRITE(*,*) ie
+                    #   WRITE(*,'(100e15.5)') x[0:4]
+                    #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+                    #   STOP "weighting 4=3=2"
+                    #
+            elif abs(x[1] - x(1)) < thr:
+                #
+                # e[3] = e[2], e[1] = e[0]
+                #
+                w2[0:2, 3] = libtetrabz_polcmplx_1221(x[3], x(2))
+                w2[0:2, 2] = w2[0:2, 3]
+                w2[0:2, 1] = libtetrabz_polcmplx_1221(x[1], x(4))
+                w2[0:2, 0] = w2[0:2, 1]
+                #
+                # IF(ANY(w2(1:2,1:4) < 0.0)):
+                #   WRITE(*,*) ie
+                #   WRITE(*,'(100e15.5)') x[0:4]
+                #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+                #   STOP "weighting 4=3 2=1"
+                #
+            else:
+                #
+                # e[3] = e[2]
+                #
+                w2[0:2, 3] = libtetrabz_polcmplx_1231(x[3], x[0], x(2))
+                w2[0:2, 2] = w2[0:2, 3]
+                w2[0:2, 1] = libtetrabz_polcmplx_1233(x[1], x[0], x(4))
+                w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[1], x(4))
+                #
+                # IF(ANY(w2(1:2,1:4) < 0.0)):
+                #   WRITE(*,*) ie
+                #   WRITE(*,'(100e15.5)') x[0:4]
+                #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+                #   STOP "weighting 4=3"
+                #
+        elif abs(x[2] - x(2)) < thr:
+            if abs(x[2] - x(1)) < thr:
+                #
+                # e[2] = e[1] = e[0]
+                #
+                w2[0:2, 3] = libtetrabz_polcmplx_1222(x[3], x(3))
+                w2[0:2, 2] = libtetrabz_polcmplx_1211(x[2], x(4))
+                w2[0:2, 1] = w2[0:2, 2]
+                w2[0:2, 0] = w2[0:2, 2]
+                #
+                # IF(ANY(w2(1:2,1:4) < 0.0)):
+                #   WRITE(*,*) ie
+                #   WRITE(*,'(100e15.5)') x[0:4]
+                #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+                #   STOP "weighting 3=2=1"
+                #
+            else:
+                #
+                # e[2] = e[1]
+                #
+                w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[0], x(3))
+                w2[0:2, 2] = libtetrabz_polcmplx_1231(x[2], x[0], x(4))
+                w2[0:2, 1] = w2[0:2, 2]
+                w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[3], x(3))
+                #
+                # IF(ANY(w2(1:2,1:4) < 0.0)):
+                #   WRITE(*,*) ie
+                #   WRITE(*,'(100e15.5)') x[0:4]
+                #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+                #   STOP "weighting 3=2"
+                #
+        elif abs(x[1] - x(1)) < thr:
+            #
+            # e[1] = e[0]
+            #
+            w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[2], x(2))
+            w2[0:2, 2] = libtetrabz_polcmplx_1233(x[2], x[3], x(2))
+            w2[0:2, 1] = libtetrabz_polcmplx_1231(x[1], x[2], x(4))
+            w2[0:2, 0] = w2[0:2, 1]
+            #
+            # IF(ANY(w2(1:2,1:4) < 0.0)):
+            #   WRITE(*,*) ie
+            #   WRITE(*,'(100e15.5)') x[0:4]
+            #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+            #   STOP "weighting 2=1"
+            #
+        else:
+            #
+            # Different each other.
+            #
+            w2[0:2, 3] = libtetrabz_polcmplx_1234(x[3], x[0], x[1], x(3))
+            w2[0:2, 2] = libtetrabz_polcmplx_1234(x[2], x[0], x[1], x(4))
+            w2[0:2, 1] = libtetrabz_polcmplx_1234(x[1], x[0], x[2], x(4))
+            w2[0:2, 0] = libtetrabz_polcmplx_1234(x[0], x[1], x[2], x(4))
+            #
+            # IF(ANY(w2(1:2,1:4) < 0.0)):
+            #   WRITE(*,*) ie
+            #   WRITE(*,'(100e15.5)') x[0:4]
+            #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+            #   STOP "weighting"
+            #
+    #
+    w1[ie, indx[0:4]] = complex(w2[0, 0:4] / e0[ie].imag, w2[1, 0:4] / (- e0[ie].imag))
+    #
+    return w1
+#
+# Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
+#  for 0<x<1, 0<y<1-x, 0<z<1-x-y
+
+
+def libtetrabz_polcmplx_1234(g1, g2, g3, g4):
+    """
+    1, Different each other
+    :param g1:
+    :param g2:
+    :param g3:
+    :param g4:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w2 = 2.0*(3.0*g2**2 - 1.0)*(math.atan(g2) - math.atan(g1)) + (g2**2 - 3.0)*g2*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = -2.0*(g2**2 - 1.0) + w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*(3.0*g3**2 - 1.0)*(math.atan(g3) - math.atan(g1)) + (g3**2 - 3.0)*g3*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = -2.0*(g3**2 - 1.0) + w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w4 = 2.0*(3.0*g4**2 - 1.0)*(math.atan(g4) - math.atan(g1)) + (g4**2 - 3.0)*g4*math.log((1.0 + g4**2)/(1.0 + g1**2))
+    w4 = -2.0*(g4**2 - 1.0) + w4/(g4 - g1)
+    w4 = w4/(g4 - g1)
+    w2 = (w2 - w3)/(g2 - g3)
+    w4 = (w4 - w3)/(g4 - g3)
+    w[0] = (w4 - w2)/(2.0*(g4 - g2))
+    #
+    # Imaginal
+    #
+    w2 = 2.0*(3.0 - g2**2)*g2*(math.atan(g2) - math.atan(g1)) + (3.0*g2**2 - 1.0)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = 4.0*g2 - w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*(3.0 - g3**2)*g3*(math.atan(g3) - math.atan(g1)) + (3.0*g3**2 - 1.0)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = 4.0*g3 - w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w4 = 2.0*(3.0 - g4**2)*g4*(math.atan(g4) - math.atan(g1)) + (3.0*g4**2 - 1.0)*math.log((1.0 + g4**2)/(1.0 + g1**2))
+    w4 = 4.0*g4 - w4/(g4 - g1)
+    w4 = w4/(g4 - g1)
+    w2 = (w2 - w3)/(g2 - g3)
+    w4 = (w4 - w3)/(g4 - g3)
+    w[1] = (w4 - w2)/(2.0*(g4 - g2))
+    #
+    return w
+
+
+def libtetrabz_polcmplx_1231(g1, g2, g3):
+    """
+    2, g4 = g1
+    :param g1:
+    :param g2:
+    :param g3:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w2 = 2.0*(-1.0 + 3.0*g2**2)*(math.atan(g2) - math.atan(g1))\
+        + g2*(-3.0 + g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = 2.0*(1.0 - g2**2) + w2/(g2 - g1)
+    w2 = -g1 + w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*(-1.0 + 3.0*g3**2)*(math.atan(g3) - math.atan(g1))\
+        + g3*(-3.0 + g3**2)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = 2.0*(1 - g3**2) + w3/(g3 - g1)
+    w3 = -g1 + w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w[0] = (w3 - w2)/(2.0*(g3 - g2))
+    #
+    # Imaginal
+    #
+    w2 = 2.0*g2*(3.0 - g2**2)*(math.atan(g2) - math.atan(g1)) + (-1.0 + 3.0*g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = 4.0*g2 - w2/(g2 - g1)
+    w2 = 1 + w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*g3*(3.0 - g3**2)*(math.atan(g3) - math.atan(g1)) + (-1.0 + 3.0*g3**2)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = 4.0*g3 - w3/(g3 - g1)
+    w3 = 1 + w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w[1] = (w3 - w2)/(2.0*(g3 - g2))
+    #
+    return w
+
+
+def libtetrabz_polcmplx_1233(g1, g2, g3):
+    """
+    3, g4 = g3
+    :param g1:
+    :param g2:
+    :param g3:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w2 = 2.0*(1.0 - 3.0*g2**2)*(math.atan(g2) - math.atan(g1)) + g2*(3.0 - g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = 2.0*(1 - g2**2) - w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*(1.0 - 3.0*g3**2)*(math.atan(g3) - math.atan(g1)) + g3*(3.0 - g3**2)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = 2.0*(1 - g3**2) - w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w2 = (w3 - w2)/(g3 - g2)
+    w3 = 4.0*(1.0 - 3.0*g1*g3)*(math.atan(g3) - math.atan(g1))\
+        + (3.0*g1 + 3.0*g3 - 3.0*g1*g3**2 + g3**3) * math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = -4.0*(1.0 - g1**2) + w3/(g3 - g1)
+    w3 = 4.0*g1 + w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w[0] = (w3 - w2)/(2.0*(g3 - g2))
+    #
+    # Imaginal
+    #
+    w2 = 2.0*g2*(3.0 - g2**2)*(math.atan(g2) - math.atan(g1)) + (-1.0 + 3.0*g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w2 = 4.0*g2 - w2/(g2 - g1)
+    w2 = w2/(g2 - g1)
+    w3 = 2.0*g3*(3.0 - g3**2)*(math.atan(g3) - math.atan(g1)) + (-1.0 + 3.0*g3**2)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = 4.0*g3 - w3/(g3 - g1)
+    w3 = w3/(g3 - g1)
+    w2 = (w3 - w2)/(g3 - g2)
+    w3 = (3.0*g1 - 3.0*g1*g3**2 + 3.0*g3 + g3**3)*(math.atan(g3) - math.atan(g1)) \
+        + (3.0*g1*g3 - 1.0)*math.log((1.0 + g3**2)/(1.0 + g1**2))
+    w3 = w3/(g3 - g1) - 4.0*g1
+    w3 = w3/(g3 - g1) - 2.0
+    w3 = (2.0*w3)/(g3 - g1)
+    w[1] = (w3 - w2)/(2.0*(g3 - g2))
+    #
+    return w
+
+
+def libtetrabz_polcmplx_1221(g1, g2):
+    """
+    4, g4 = g1 and g3 = g2
+    :param g1:
+    :param g2:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w[0] = -2.0*(-1.0 + 2.0*g1*g2 + g2**2)*(math.atan(g2) - math.atan(g1)) \
+        + (g1 + 2.0*g2 - g1*g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w[0] = 2.0*(-1.0 + g1**2) + w[0]/(g2 - g1)
+    w[0] = 3.0*g1 + w[0]/(g2 - g1)
+    w[0] = 2.0 + (3.0*w[0])/(g2 - g1)
+    w[0] = w[0]/(2.0*(g2 - g1))
+    #
+    # Imaginal
+    #
+    w[1] = 2.0*(g1 + 2.0*g2 - g1*g2**2)*(math.atan(g2) - math.atan(g1)) \
+        + (-1.0 + 2.0*g1*g2 + g2**2)*math.log((1 + g2**2)/(1 + g1**2))
+    w[1] = -4.0*g1 + w[1]/(g2 - g1)
+    w[1] = -3.0 + w[1]/(g2 - g1)
+    w[1] = (3.0*w[1])/(2.0*(g2 - g1)**2)
+    #
+    return w
+
+
+def libtetrabz_polcmplx_1222(g1, g2):
+    """
+    5, g4 = g3 = g2
+    :param g1:
+    :param g2:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w[0] = 2.0*(-1.0 + g1**2 + 2.0*g1*g2)*(math.atan(g2) - math.atan(g1)) \
+        + (-2.0*g1 - g2 + g1**2*g2) * math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w[0] = 2.0*(1.0 - g1**2) + w[0]/(g2 - g1)
+    w[0] = g1 - w[0]/(g2 - g1)
+    w[0] = 1.0 - (3.0*w[0])/(g2 - g1)
+    w[0] = w[0]/(2.0*(g2 - g1))
+    #
+    # Imaginal
+    #
+    w[1] = 2.0*(-2.0*g1 - g2 + g1**2*g2)*(math.atan(g2) - math.atan(g1)) \
+        + (1.0 - g1**2 - 2.0*g1*g2) * math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w[1] = 4.0*g1 + w[1]/(g2 - g1)
+    w[1] = 1.0 + w[1]/(g2 - g1)
+    w[1] = (3.0*w[1])/(2.0*(g2 - g1)**2)
+    #
+    return w
+
+
+def libtetrabz_polcmplx_1211(g1, g2):
+    """
+    6, g4 = g3 = g1
+    :param g1:
+    :param g2:
+    :return:
+    """
+    w = numpy.empty(2, dtype=numpy.float_)
+    #
+    # Real
+    #
+    w[0] = 2.0*(3.0*g2**2 - 1.0)*(math.atan(g2) - math.atan(g1)) \
+        + g2*(g2**2 - 3.0)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w[0] = 2.0*(1.0 - g1**2) + w[0]/(g2 - g1)
+    w[0] = -5.0*g1 + w[0]/(g2 - g1)
+    w[0] = -11.0 + (3.0*w[0])/(g2 - g1)
+    w[0] = w[0]/(6.0*(g2 - g1))
+    #
+    # Imaginal
+    #
+    w[1] = 2.0*g2*(-3.0 + g2**2)*(math.atan(g2) - math.atan(g1)) \
+        + (1.0 - 3.0*g2**2)*math.log((1.0 + g2**2)/(1.0 + g1**2))
+    w[1] = 4.0*g2 + w[1]/(g2 - g1)
+    w[1] = 1.0 + w[1]/(g2 - g1)
+    w[1] = w[1]/(2.0*(g2 - g1)**2)
+    #
+    return w
similarity index 99%
rename from src/libtetrabz/libtetrabz_polstat_mod.py
rename to src/libtetrabz/libtetrabz_polstat.py
index 5b9fb95..be04f43 100644 (file)
@@ -25,7 +25,7 @@ import numpy
 import libtetrabz_common
 
 
-def libtetrabz_polstat(bvec, eig1, eig2):
+def polstat(bvec, eig1, eig2):
     """
     Compute Static polarization function
     :param bvec:
@@ -134,7 +134,7 @@ def libtetrabz_polstat(bvec, eig1, eig2):
                 ei2 = ei1[0:4,   ib]
                 ej2 = ej1[0:4, 0:nb]
                 w2 = libtetrabz_polstat2(ei2, ej2)
-                w1[0:nb, 0:4] += w2[1:nb, 0:4]
+                w1[0:nb, 0:4] += w2[0:nb, 0:4]
                 #
             else:
                 continue
@@ -228,7 +228,7 @@ def libtetrabz_polstat2(ei1, ej1):
             #
             de = ej1[0:4, ib] - ei1[0:4]
             w2 = libtetrabz_polstat3(de)
-            w1[ib, indx[0:4]] += w2
+            w1[ib, 0:4] += w2
     #
     return w1