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