1#!/usr/local/bin/python3.8
2# irrbz.py
3#
4# Usage: irrbz.py irrbz prefix
5#
6# Extracts irreducible wedge from siesta output files prefix.KP
7# and prefix.EIG and generates prefix.KP.irrbz and prefix.EIG.irrbz.
8# Reads in lattice vectors from prefix.XV and k-points from irrbz
9# generated by kgrid.x.
10#
11# Georgy Samsonidze, UCB, December 2010
12
13eps = 1.0e-5
14
15def main(argv = None):
16   if argv is None:
17      argv = sys.argv
18   argc = len(argv)
19   self = "irrbz.py"
20   if argv[0][-len(self):] != self:
21      sys.exit("\n   Rename script to %s\n" % self)
22   if argc != 3:
23      sys.exit("\n   Usage: %s irrbz prefix\n" % self)
24   irrbz = argv[1]
25   prefix = argv[2]
26   kpin = prefix + '.KP'
27   kpout = prefix + '.KP.irrbz'
28   eigin = prefix + '.EIG'
29   eigout = prefix + '.EIG.irrbz'
30
31   h = open(irrbz, 'r')
32   s = h.readline()
33   if s[0:16] != 'K_POINTS crystal':
34      sys.exit("\n   Error: irrbz file %s\n" % irrbz)
35   s = h.readline()
36   nirr = int(s)
37   kirr = []
38   for i in xrange(nirr):
39      s = h.readline()
40      t = s.split()
41      k1 = float(t[0])
42      k2 = float(t[1])
43      k3 = float(t[2])
44      kirr.append([k1,k2,k3])
45   h.close()
46
47   h = open(prefix + '.XV', 'r')
48   avec = []
49   for i in xrange(3):
50      s = h.readline()
51      t = s.split()
52      ax = float(t[0])
53      ay = float(t[1])
54      az = float(t[2])
55      avec.append([ax,ay,az])
56   h.close()
57   bvec = [[], [], []]
58   celvol = 0.0
59   celvol += avec[0][0] * (avec[1][1] * avec[2][2] - avec[2][1] * avec[1][2])
60   celvol -= avec[0][1] * (avec[1][0] * avec[2][2] - avec[2][0] * avec[1][2])
61   celvol += avec[0][2] * (avec[1][0] * avec[2][1] - avec[2][0] * avec[1][1])
62   bvec[0].append(2.0 * math.pi * (avec[1][1] * avec[2][2] - avec[1][2] * avec[2][1]) / celvol)
63   bvec[0].append(2.0 * math.pi * (avec[1][2] * avec[2][0] - avec[1][0] * avec[2][2]) / celvol)
64   bvec[0].append(2.0 * math.pi * (avec[1][0] * avec[2][1] - avec[1][1] * avec[2][0]) / celvol)
65   bvec[1].append(2.0 * math.pi * (avec[2][1] * avec[0][2] - avec[2][2] * avec[0][1]) / celvol)
66   bvec[1].append(2.0 * math.pi * (avec[2][2] * avec[0][0] - avec[2][0] * avec[0][2]) / celvol)
67   bvec[1].append(2.0 * math.pi * (avec[2][0] * avec[0][1] - avec[2][1] * avec[0][0]) / celvol)
68   bvec[2].append(2.0 * math.pi * (avec[0][1] * avec[1][2] - avec[0][2] * avec[1][1]) / celvol)
69   bvec[2].append(2.0 * math.pi * (avec[0][2] * avec[1][0] - avec[0][0] * avec[1][2]) / celvol)
70   bvec[2].append(2.0 * math.pi * (avec[0][0] * avec[1][1] - avec[0][1] * avec[1][0]) / celvol)
71
72   h = open(kpin, 'r')
73   s = h.readline()
74   nful = int(s)
75   kful = []
76   for i in xrange(nful):
77      s = h.readline()
78      t = s.split()
79      j = int(t[0])
80      kx = float(t[1])
81      ky = float(t[2])
82      kz = float(t[3])
83      wk = float(t[4])
84      if j != i + 1:
85         sys.exit("\n   Error: KP file %s\n" % kpin)
86      kful.append([kx,ky,kz,wk])
87   h.close()
88
89   kidx = []
90   for i in xrange(nirr):
91      for j in xrange(nful):
92         kx = kful[j][0]
93         ky = kful[j][1]
94         kz = kful[j][2]
95         k1 = (kx * avec[0][0] + ky * avec[0][1] + kz * avec[0][2]) / (2.0 * math.pi)
96         k2 = (kx * avec[1][0] + ky * avec[1][1] + kz * avec[1][2]) / (2.0 * math.pi)
97         k3 = (kx * avec[2][0] + ky * avec[2][1] + kz * avec[2][2]) / (2.0 * math.pi)
98         d1 = kirr[i][0] - k1
99         d2 = kirr[i][1] - k2
100         d3 = kirr[i][2] - k3
101         l1 = abs(d1 - round(d1)) + abs(d2 - round(d2)) + abs(d3 - round(d3))
102         d1 = kirr[i][0] + k1
103         d2 = kirr[i][1] + k2
104         d3 = kirr[i][2] + k3
105         l2 = abs(d1 - round(d1)) + abs(d2 - round(d2)) + abs(d3 - round(d3))
106         if l1 < eps or l2 < eps:
107            kidx.append(j)
108            break
109   if len(kidx) != nirr:
110      sys.exit("\n   Error: cannot match k-points\n")
111
112   kout = []
113   for i in xrange(nirr):
114      k1 = kirr[i][0]
115      k2 = kirr[i][1]
116      k3 = kirr[i][2]
117      kx = (k1 * bvec[0][0] + k2 * bvec[1][0] + k3 * bvec[2][0])
118      ky = (k1 * bvec[0][1] + k2 * bvec[1][1] + k3 * bvec[2][1])
119      kz = (k1 * bvec[0][2] + k2 * bvec[1][2] + k3 * bvec[2][2])
120      wk = kful[kidx[i]][3]
121      kout.append([kx,ky,kz,wk])
122
123   h = open(kpout, 'w')
124   s = '%6i\n' % nirr
125   h.write(s)
126   for i in xrange(nirr):
127      s = '%6i%12.6f%12.6f%12.6f%15.6f\n' % (i + 1, kout[i][0], kout[i][1], kout[i][2], kout[i][3])
128      h.write(s)
129   h.close()
130
131   h = open(eigin, 'r')
132   s = h.readline()
133   efermi = float(s)
134   s = h.readline()
135   t = s.split()
136   nb = int(t[0])
137   ns = int(t[1])
138   nk = int(t[2])
139   if nk != nful:
140      sys.exit("\n   Error: EIG file %s\n" % eigin)
141   nline = (nb * ns) / 10
142   if 10 * nline < nb * ns:
143      nline += 1
144   en = []
145   for i in xrange(nful):
146      en.append([])
147      for j in xrange(nline):
148         s = h.readline()
149         if j == 0:
150            t = s.split()
151            if int(t[0]) != i + 1:
152               sys.exit("\n   Error: EIG file %s\n" % eigin)
153         en[i].append(s)
154   h.close()
155
156   h = open(eigout, 'w')
157   s = '%14.4f\n' % efermi
158   h.write(s)
159   s = '%6i%6i%6i\n' % (nb, ns, nirr)
160   h.write(s)
161   for i in xrange(nirr):
162      for j in xrange(nline):
163         if j == 0:
164            s = '%5i' % (i + 1)
165            s += en[kidx[i]][j][5:]
166            h.write(s)
167         else:
168            h.write(en[kidx[i]][j])
169   h.close()
170
171   return 0
172
173if __name__ == "__main__":
174   import sys
175   import math
176   sys.exit(main())
177