get rid of np.matrix

nanonet/negf/greens_functions.py

@@ -73,8 +73,7 @@ def surface_greens_function_poles(h_list):
     eigenvects = eigenvects[:, ind]
 
     eigenvects = eigenvects[matix_size:, :]
-    eigenvals = np.matrix(np.diag(eigenvals))
-    eigenvects = np.matrix(eigenvects)
+    eigenvals = np.diag(eigenvals)
 
     norms = linalg.norm(eigenvects, axis=0)
     norms = np.array([1e30 if np.abs(norm) < 0.000001 else norm for norm in norms])
@@ -97,13 +96,13 @@ def group_velocity(eigenvector, eigenvalue, h_r):
                               an eigenvector and an eigenvalue
     """
 
-    return np.imag(eigenvector.H * h_r * eigenvalue * eigenvector)
+    return np.imag(np.dot(np.dot(np.dot(eigenvector.conj().T, h_r), eigenvalue), eigenvector))
 
 
 def iterate_gf(E, h_0, h_l, h_r, gf, num_iter):
 
     for _ in range(num_iter):
-        gf = h_r * np.linalg.pinv(E * np.identity(h_0.shape[0]) - h_0 - gf) * h_l
+        gf = h_r.dot(np.linalg.pinv(E * np.identity(h_0.shape[0]) - h_0 - gf)).dot(h_l)
 
     return gf
 
@@ -129,10 +128,10 @@ def surface_greens_function(E, h_l, h_0, h_r, iterate=True, damp=0.0001j):
     vals, vects = surface_greens_function_poles(h_list)
     vals = np.diag(vals)
 
-    u_right = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    u_left = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    lambda_right = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    lambda_left = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
+    u_right = np.zeros(h_0.shape, dtype=np.complex)
+    u_left = np.zeros(h_0.shape, dtype=np.complex)
+    lambda_right = np.zeros(h_0.shape, dtype=np.complex)
+    lambda_left = np.zeros(h_0.shape, dtype=np.complex)
 
     alpha = 0.01
 
@@ -171,8 +170,8 @@ def surface_greens_function(E, h_l, h_0, h_r, iterate=True, damp=0.0001j):
                 lambda_left[j, j] = vals[-j + 2*h_0.shape[0]-1]
                 u_left[:, j] = vects[:, -j + 2*h_0.shape[0]-1]
 
-    sgf_l = h_r * u_right * lambda_right * np.linalg.pinv(u_right)
-    sgf_r = h_l * u_left * lambda_right * np.linalg.pinv(u_left)
+    sgf_l = h_r.dot(u_right).dot(lambda_right).dot(np.linalg.pinv(u_right))
+    sgf_r = h_l.dot(u_left).dot(lambda_right).dot(np.linalg.pinv(u_left))
 
     if iterate:
         return iterate_gf(E, h_0, h_l, h_r, sgf_l, 2), iterate_gf(E, h_0, h_r, h_l, sgf_r, 2)

nanonet/negf/hamiltonian_chain.py

@@ -40,11 +40,11 @@ def sgf(self):
 
         for jjj in range(len(self.h_0)):
             if jjj == 0:
-                sgf[jjj] = -2.0 * np.matrix(np.imag(self.sgf_r) * fd(self.energy, self.ef1, self.tempr))
+                sgf[jjj] = -2.0 * np.imag(self.sgf_r) * fd(self.energy, self.ef1, self.tempr)
             elif jjj == len(self.h_0) - 1:
-                sgf[jjj] = -2.0 * np.matrix(np.imag(self.sgf_l) * fd(self.energy, self.ef2, self.tempr))
+                sgf[jjj] = -2.0 * np.imag(self.sgf_l) * fd(self.energy, self.ef2, self.tempr)
             else:
-                sgf[jjj] = np.matrix(np.zeros(self.h_0[jjj].shape))
+                sgf[jjj] = np.zeros(self.h_0[jjj].shape)
 
         return sgf
 

nanonet/negf/recursive_greens_functions.py

@@ -77,20 +77,20 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
     gr_left[0] = mat_left_div(-mat_d_list[0], np.eye(mat_shapes[0][0]))  # Initialising the retarded left connected.
 
     for q in range(num_of_matrices - 1):  # Recursive algorithm (B2)
-        gr_left[q + 1] = mat_left_div((-mat_d_list[q + 1] - mat_l_list[q] * gr_left[q] * mat_u_list[q]),
+        gr_left[q + 1] = mat_left_div((-mat_d_list[q + 1] - mat_l_list[q].dot(gr_left[q]).dot(mat_u_list[q])),
                                       np.eye(mat_shapes[q + 1][0]))      # The left connected recursion.
     # -------------------------------------------------------------------
 
     grl = [None for _ in range(num_of_matrices-1)]
     gru = [None for _ in range(num_of_matrices-1)]
-    grd = copy.copy(gr_left)                                       # Our glorious benefactor.
+    grd = copy.copy(gr_left)                                             # Our glorious benefactor.
     g_trans = copy.copy(gr_left[len(gr_left)-1])
 
-    for q in range(num_of_matrices - 2, -1, -1):                   # Recursive algorithm
-        grl[q] = grd[q + 1] * mat_l_list[q] * gr_left[q]           # (B5) We get the off-diagonal blocks for free.
-        gru[q] = gr_left[q] * mat_u_list[q] * grd[q + 1]           # (B6) because we need .Tthem.T for the next calc:
-        grd[q] = gr_left[q] + gr_left[q] * mat_u_list[q] * grl[q]  # (B4) I suppose I could also use the lower.
-        g_trans = gr_left[q] * mat_u_list[q] * g_trans
+    for q in range(num_of_matrices - 2, -1, -1):                         # Recursive algorithm
+        grl[q] = grd[q + 1].dot(mat_l_list[q]).dot(gr_left[q])           # (B5) We get the off-diagonal blocks for free.
+        gru[q] = gr_left[q].dot(mat_u_list[q]).dot(grd[q + 1])           # (B6) because we need .Tthem.T for the next calc:
+        grd[q] = gr_left[q] + gr_left[q].dot(mat_u_list[q]).dot(grl[q])  # (B4) I suppose I could also use the lower.
+        g_trans = gr_left[q].dot(mat_u_list[q]).dot(g_trans)
 
     # -------------------------------------------------------------------
     # ------ compute the electron correlation function if needed --------
@@ -99,12 +99,12 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
     if isinstance(s_in, list):
 
         gin_left = [None for _ in range(num_of_matrices)]
-        gin_left[0] = gr_left[0] * s_in[0] * np.conj(gr_left[0])
+        gin_left[0] = gr_left[0].dot(s_in[0]).dot(np.conj(gr_left[0]))
 
         for q in range(num_of_matrices - 1):
-            sla2 = mat_l_list[q] * gin_left[q] * np.conj(mat_u_list[q])
+            sla2 = mat_l_list[q].dot(gin_left[q]).dot(np.conj(mat_u_list[q]))
             prom = s_in[q + 1] + sla2
-            gin_left[q + 1] = np.real(gr_left[q + 1] * prom * np.conj(gr_left[q + 1]))
+            gin_left[q + 1] = np.real(gr_left[q + 1].dot(prom).dot(np.conj(gr_left[q + 1])))
 
         # ---------------------------------------------------------------
 
@@ -113,10 +113,11 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
         gnd = copy.copy(gin_left)
 
         for q in range(num_of_matrices - 2, -1, -1):               # Recursive algorithm
-            gnl[q] = grd[q + 1] * mat_l_list[q] * gin_left[q] + gnd[q + 1] * np.conj(mat_l_list[q]) * np.conj(gr_left[q])
+            gnl[q] = grd[q + 1].dot(mat_l_list[q]).dot(gin_left[q]) +\
+                     gnd[q + 1].dot(np.conj(mat_l_list[q])).dot(np.conj(gr_left[q]))
             gnd[q] = np.real(gin_left[q] +
-                             gr_left[q] * mat_u_list[q] * gnd[q + 1] * np.conj(mat_l_list[q]) * np.conj(gr_left[q]) +
-                             (gin_left[q] * np.conj(mat_u_list[q]) * np.conj(grl[q]) + gru[q] * mat_l_list[q] * gin_left[q]))
+                             gr_left[q].dot(mat_u_list[q]).dot(gnd[q + 1]).dot(np.conj(mat_l_list[q])).dot(np.conj(gr_left[q])) +
+                             (gin_left[q].dot(np.conj(mat_u_list[q])).dot(np.conj(grl[q])) + gru[q].dot(mat_l_list[q]).dot(gin_left[q])))
 
             gnu[q] = np.conj(gnl[q])
 
@@ -126,12 +127,12 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
     if isinstance(s_out, list):
 
         gip_left = [None for _ in range(num_of_matrices)]
-        gip_left[0] = gr_left[0] * s_out[0] * np.conj(gr_left[0])
+        gip_left[0] = gr_left[0].dot(s_out[0]).dot(np.conj(gr_left[0]))
 
         for q in range(num_of_matrices - 1):
-            sla2 = mat_l_list[q] * gip_left[q] * np.conj(mat_u_list[q])
+            sla2 = mat_l_list[q].dot(gip_left[q]).dot(np.conj(mat_u_list[q]))
             prom = s_out[q + 1] + sla2
-            gip_left[q + 1] = np.real(gr_left[q + 1] * prom * np.conj(gr_left[q + 1]))
+            gip_left[q + 1] = np.real(gr_left[q + 1].dot(prom).dot(np.conj(gr_left[q + 1])))
 
         # ---------------------------------------------------------------
 
@@ -140,12 +141,12 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
         gpd = copy.copy(gip_left)
 
         for q in range(num_of_matrices - 2, -1, -1):               # Recursive algorithm
-            gpl[q] = grd[q + 1] * mat_l_list[q] * gip_left[q] + gpd[q + 1] * np.conj(mat_l_list[q]) * np.conj(gr_left[q])
+            gpl[q] = grd[q + 1].dot(mat_l_list[q]).dot(gip_left[q]) + gpd[q + 1].dot(np.conj(mat_l_list[q])).dot(np.conj(gr_left[q]))
             gpd[q] = np.real(gip_left[q] +
-                             gr_left[q] * mat_u_list[q] * gpd[q + 1] * np.conj(mat_l_list[q]) * np.conj(gr_left[q]) +
-                             (gip_left[q] * np.conj(mat_u_list[q]) * np.conj(grl[q]) + gru[q] * mat_l_list[q] * gip_left[q]))
+                             gr_left[q].dot(mat_u_list[q]).dot(gpd[q + 1]).dot(np.conj(mat_l_list[q])).dot(np.conj(gr_left[q])) +
+                             gip_left[q].dot(np.conj(mat_u_list[q])).dot(np.conj(grl[q])) + gru[q].dot(mat_l_list[q]).dot(gip_left[q]))
 
-            gpu[0] = gpl[0].H
+            gpu[0] = gpl[0].conj().T
 
     # -------------------------------------------------------------------
     # -- remove energy from the main diagonal of th Hamiltonian matrix --

nanonet/negf/surf_greens_function.py

@@ -44,13 +44,13 @@ def sort_eigs(alpha, betha, Z, h_l, h_r, flag):
     inds_deg = [i for i in inds_deg if len(i) > 1]
 
     for item in inds_deg:
-        phi = np.matrix(Z[h_r.shape[0]:, item])
-        oper = 1j * np.matrix(h_r * eigenvals[item[0]] - h_l * np.conj(eigenvals[item[0]]))
-        _, vects = np.linalg.eig(phi.H * oper * phi)
+        phi = Z[h_r.shape[0]:, item]
+        oper = 1j * np.dot(h_r, eigenvals[item[0]]) - np.dot(h_l, eigenvals[item[0]].conj())
+        _, vects = np.linalg.eig(np.dot(np.dot(phi.conj().T, oper), phi))
         print(vects)
-        Z[h_r.shape[0]:, item] = phi * np.matrix(vects)
+        Z[h_r.shape[0]:, item] = np.dot(phi, vects)
 
-    phi = np.matrix(Z)
+    phi = Z
 
     for j, item in enumerate(zip(alpha, betha)):
 
@@ -64,7 +64,7 @@ def sort_eigs(alpha, betha, Z, h_l, h_r, flag):
         elif np.abs(eigenval) < 1.0 - margin:
             ans.append(True)
         else:
-            gv = np.imag(2 * phi[h_r.shape[0]:, j].H * np.matrix(h_r) * phi[h_r.shape[0]:, j])
+            gv = np.imag(2 * np.dot(np.dot(phi[h_r.shape[0]:, j].conj().T, h_r), phi[h_r.shape[0]:, j]))
             # gv = group_velocity(phi[h_r.shape[0]:, j], eigenval, h_r)
 
             if flag:
@@ -146,8 +146,7 @@ def surface_greens_function_poles(E, h_l, h_0, h_r):
     eigenvects = eigenvects[:, ind]
 
     eigenvects = eigenvects[h_0.shape[0]:, :]
-    eigenvals = np.matrix(np.diag(eigenvals))
-    eigenvects = np.matrix(eigenvects)
+    eigenvals = np.diag(eigenvals)
 
     norms = linalg.norm(eigenvects, axis=0)
     norms = np.array([1e30 if np.abs(norm) < 0.000001 else norm for norm in norms])
@@ -193,11 +192,8 @@ def surface_greens_function_poles_Shur(E, h_l, h_0, h_r):
                                                                      output='complex',
                                                                      sort=sort1)
 
-    eigv_left = np.matrix(eigv_left)
-    eigv_left1 = np.matrix(eigv_left1)
-
-    return h_r * eigv_left[h_0.shape[0]:, :h_0.shape[0]] * np.linalg.pinv(eigv_left[:h_0.shape[0], :h_0.shape[0]]), \
-           h_l * eigv_left1[h_0.shape[0]:, :h_0.shape[0]] * np.linalg.pinv(eigv_left1[:h_0.shape[0], :h_0.shape[0]])
+    return h_r.dot(eigv_left[h_0.shape[0]:, :h_0.shape[0]]).dot(np.linalg.pinv(eigv_left[:h_0.shape[0], :h_0.shape[0]])), \
+           h_l.dot(eigv_left1[h_0.shape[0]:, :h_0.shape[0]]).dot(np.linalg.pinv(eigv_left1[:h_0.shape[0], :h_0.shape[0]]))
 
 
 def group_velocity(eigenvector, eigenvalue, h_r):
@@ -210,7 +206,7 @@ def group_velocity(eigenvector, eigenvalue, h_r):
     :return:
     """
 
-    return np.imag(eigenvector.H * h_r * eigenvalue * eigenvector)
+    return np.imag(np.dot(np.dot(np.dot(eigenvector.conj().T, h_r), eigenvalue), eigenvector))
 
 
 def surface_greens_function(E, h_l, h_0, h_r):
@@ -229,10 +225,10 @@ def surface_greens_function(E, h_l, h_0, h_r):
     vals, vects = surface_greens_function_poles(E, h_l, h_0, h_r)
     vals = np.diag(vals)
 
-    u_right = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    u_left = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    lambda_right = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
-    lambda_left = np.matrix(np.zeros(h_0.shape, dtype=np.complex))
+    u_right = np.zeros(h_0.shape, dtype=np.complex)
+    u_left = np.zeros(h_0.shape, dtype=np.complex)
+    lambda_right = np.zeros(h_0.shape, dtype=np.complex)
+    lambda_left = np.zeros(h_0.shape, dtype=np.complex)
 
     alpha = 0.001
 
@@ -279,8 +275,8 @@ def surface_greens_function(E, h_l, h_0, h_r):
             # lambda_left[j, j] = vals[-j + 2 * h_0.shape[0] - 1]
             # u_left[:, j] = vects[:, -j + 2 * h_0.shape[0] - 1]
 
-    sgf_l = h_r * u_right * lambda_right * np.linalg.pinv(u_right)
-    sgf_r = h_l * u_left * lambda_right * np.linalg.pinv(u_left)
+    sgf_l = h_r.dot(u_right).dot(lambda_right).dot(np.linalg.pinv(u_right))
+    sgf_r = h_l.dot(u_left).dot(lambda_right).dot(np.linalg.pinv(u_left))
 
     # sgf_l = u_right[h_0.shape[0]:, :] * np.linalg.pinv(u_right[:h_0.shape[0], :])
     # sgf_r = h_l * u_left * lambda_right * np.linalg.pinv(u_left)
@@ -298,7 +294,7 @@ def surface_greens_function(E, h_l, h_0, h_r):
 def main():
     import sys
     sys.path.insert(0, '/home/mk/TB_project/tb')
-    import tb
+    import nanonet.tb as tb
 
     a = tb.Atom('A')
     a.add_orbital('s', -0.7)
@@ -364,18 +360,18 @@ def main():
     tr = np.zeros((energy.shape[0]), dtype=np.complex)
 
     for j, E in enumerate(energy):
-        gf0 = np.matrix(gf[j, :, :])
-        gamma_l = 1j * (np.matrix(sgf_l[j, :, :]) - np.matrix(sgf_l[j, :, :]).H)
-        gamma_r = 1j * (np.matrix(sgf_r[j, :, :]) - np.matrix(sgf_r[j, :, :]).H)
-        tr[j] = np.real(np.trace(gamma_l * gf0 * gamma_r * gf0.H))
-        dos[j] = np.real(np.trace(1j * (gf0 - gf0.H)))
+        gf0 = gf[j, :, :]
+        gamma_l = 1j * (sgf_l[j, :, :] - sgf_l[j, :, :].conj().T)
+        gamma_r = 1j * (sgf_r[j, :, :] - sgf_r[j, :, :].conj().T)
+        tr[j] = np.real(np.trace(gamma_l.dot(gf0).dot(gamma_r).dot(gf0.conj().T)))
+        dos[j] = np.real(np.trace(1j * (gf0 - gf0.conj().T)))
     print(sgf_l.shape)
 
 
 def main1():
     import sys
     sys.path.insert(0, '/home/mk/TB_project/tb')
-    import tb
+    import nanonet.tb as tb
 
     tb.Atom.orbital_sets = {'Si': 'SiliconSP3D5S', 'H': 'HydrogenS'}
     h = tb.Hamiltonian(xyz='/home/mk/NEGF_project/SiNW.xyz', nn_distance=2.4)
@@ -434,11 +430,11 @@ def main1():
     dos = np.zeros((energy.shape[0]), dtype=np.complex)
 
     for j, E in enumerate(energy):
-        gf0 = np.matrix(gf[j, :, :])
-        gamma_l = 1j * (np.matrix(sgf_l[j, :, :]) - np.matrix(sgf_l[j, :, :]).H)
-        gamma_r = 1j * (np.matrix(sgf_r[j, :, :]) - np.matrix(sgf_r[j, :, :]).H)
-        dos[j] = np.real(np.trace(1j * (gf0 - gf0.H)))
-        tr[j] = np.real(np.trace(gamma_l * gf0 * gamma_r * gf0.H))
+        gf0 = gf[j, :, :]
+        gamma_l = 1j * (sgf_l[j, :, :] - sgf_l[j, :, :].conj().T)
+        gamma_r = 1j * (sgf_r[j, :, :] - sgf_r[j, :, :].conj().T)
+        dos[j] = np.real(np.trace(1j * (gf0 - gf0.conj().T)))
+        tr[j] = np.real(np.trace(gamma_l.dot(gf0).dot(gamma_r).dot(gf0.conj().T)))
 
     ax = plt.axes()
     ax.set_xlabel('Energy (eV)')
@@ -477,7 +473,7 @@ def regularize_gf(gf):
 def inverse_bs_problem():
     import sys
     sys.path.insert(0, '/home/mk/TB_project/tb')
-    import tb
+    import nanonet.tb as tb
 
     # a = tb.Atom('A')
     # a.add_orbital('s', -0.7)
@@ -553,7 +549,7 @@ def inverse_bs_problem():
 def main2():
     import sys
     sys.path.insert(0, '/home/mk/TB_project/tb')
-    import tb
+    import nanonet.tb as tb
 
     a = tb.Atom('A')
     a.add_orbital('s', -0.7)
@@ -598,18 +594,18 @@ def main2():
     tr = np.zeros((energy.shape[0]), dtype=np.complex)
 
     for j, E in enumerate(energy):
-        gf0 = np.matrix(gf[j, :, :])
-        gamma_l = 1j * (np.matrix(sgf_l[j, :, :]) - np.matrix(sgf_l[j, :, :]).H)
-        gamma_r = 1j * (np.matrix(sgf_r[j, :, :]) - np.matrix(sgf_r[j, :, :]).H)
-        tr[j] = np.real(np.trace(gamma_l * gf0 * gamma_r * gf0.H))
-        dos[j] = np.real(np.trace(1j * (gf0 - gf0.H)))
+        gf0 = gf[j, :, :]
+        gamma_l = 1j * (sgf_l[j, :, :] - sgf_l[j, :, :].conj().T)
+        gamma_r = 1j * (sgf_r[j, :, :] - sgf_r[j, :, :].conj().T)
+        tr[j] = np.real(np.trace(gamma_l.dot(gf0).dot(gamma_r).dot(gf0.conj().T)))
+        dos[j] = np.real(np.trace(1j * (gf0 - gf0.conj().T)))
     print(sgf_l.shape)
 
 
 def main3():
     import sys
     sys.path.insert(0, '/home/mk/TB_project/tb')
-    import tb
+    import nanonet.tb as tb
 
     a = tb.Atom('A')
     a.add_orbital('s', -0.7)
@@ -654,11 +650,11 @@ def main3():
     tr = np.zeros((energy.shape[0]), dtype=np.complex)
 
     for j, E in enumerate(energy):
-        gf0 = np.matrix(gf[j, :, :])
-        gamma_l = 1j * (np.matrix(sgf_l[j, :, :]) - np.matrix(sgf_l[j, :, :]).H)
-        gamma_r = 1j * (np.matrix(sgf_r[j, :, :]) - np.matrix(sgf_r[j, :, :]).H)
-        tr[j] = np.real(np.trace(gamma_l * gf0 * gamma_r * gf0.H))
-        dos[j] = np.real(np.trace(1j * (gf0 - gf0.H)))
+        gf0 = gf[j, :, :]
+        gamma_l = 1j * (sgf_l[j, :, :] - sgf_l[j, :, :].conj().T)
+        gamma_r = 1j * (sgf_r[j, :, :] - sgf_r[j, :, :].conj().T)
+        tr[j] = np.real(np.trace(gamma_l.dot(gf0).dot(gamma_r).dot(gf0.conj().T)))
+        dos[j] = np.real(np.trace(1j * (gf0 - gf0.conj().T)))
 
     print(sgf_l.shape)
 

nanonet/tb/constants.py

@@ -31,28 +31,28 @@
 b_bi = a_bi / c_bi
 gamma_bi = 0.2303
 
-recip_lat_basis = np.matrix([[-1, -np.sqrt(3.0)/3.0, b_bi],
+recip_lat_basis = np.array([[-1, -np.sqrt(3.0)/3.0, b_bi],
                              [1, -np.sqrt(3.0)/3.0, b_bi],
                              [0, 2*np.sqrt(3.0)/3.0, b_bi]]) * 2 * PI / a_bi
 
 SPECIAL_K_POINTS_BI = {
-    'GAMMA': (recip_lat_basis * np.matrix([0, 0, 0]).T).tolist(),
-    'X': (recip_lat_basis * np.matrix([0.5, 0.5, 0.0]).T).tolist(),
-    'T': (recip_lat_basis * np.matrix([0.5, 0.5, 0.5]).T).tolist(),
-    'L': (recip_lat_basis * np.matrix([0.0, 0.5, 0.0]).T).tolist(),
-    'W': (recip_lat_basis * np.matrix([gamma_bi, 1.0-gamma_bi, 0.5]).T).tolist(),
-    'K': (recip_lat_basis * np.matrix([1.0/4+0.5*gamma_bi, 3.0/4-0.5*gamma_bi, 0.0]).T).tolist(),
-    'U': (recip_lat_basis * np.matrix([0.5*gamma_bi+0.25, 1.0-gamma_bi, 0.5*gamma_bi+0.25]).T).tolist()
+    'GAMMA': (recip_lat_basis * np.array([0, 0, 0]).T).tolist(),
+    'X': (recip_lat_basis * np.array([0.5, 0.5, 0.0]).T).tolist(),
+    'T': (recip_lat_basis * np.array([0.5, 0.5, 0.5]).T).tolist(),
+    'L': (recip_lat_basis * np.array([0.0, 0.5, 0.0]).T).tolist(),
+    'W': (recip_lat_basis * np.array([gamma_bi, 1.0-gamma_bi, 0.5]).T).tolist(),
+    'K': (recip_lat_basis * np.array([1.0/4+0.5*gamma_bi, 3.0/4-0.5*gamma_bi, 0.0]).T).tolist(),
+    'U': (recip_lat_basis * np.array([0.5*gamma_bi+0.25, 1.0-gamma_bi, 0.5*gamma_bi+0.25]).T).tolist()
 }
 
 # aaa = 3.289
 #
 # SPECIAL_K_POINTS_BI = {
 #     'GAMMA': [0, 0, 0],
 #     'X': [0.5 * PI/aaa, 0.5 * PI/aaa, 0.5 * PI/aaa],
-#     'T': (recip_lat_basis * np.matrix([0.5, 0.5, 0.5]).T).tolist(),
-#     'L': (recip_lat_basis * np.matrix([0.0, 0.5, 0.0]).T).tolist(),
-#     'W': (recip_lat_basis * np.matrix([gamma_bi, 1.0-gamma_bi, 0.5]).T).tolist(),
-#     'K': (recip_lat_basis * np.matrix([1.0/4+0.5*gamma_bi, 3.0/4-0.5*gamma_bi, 0.0]).T).tolist(),
-#     'U': (recip_lat_basis * np.matrix([0.5*gamma_bi+0.25, 1.0-gamma_bi, 0.5*gamma_bi+0.25]).T).tolist()
+#     'T': (recip_lat_basis * np.array([0.5, 0.5, 0.5]).T).tolist(),
+#     'L': (recip_lat_basis * np.array([0.0, 0.5, 0.0]).T).tolist(),
+#     'W': (recip_lat_basis * np.array([gamma_bi, 1.0-gamma_bi, 0.5]).T).tolist(),
+#     'K': (recip_lat_basis * np.array([1.0/4+0.5*gamma_bi, 3.0/4-0.5*gamma_bi, 0.0]).T).tolist(),
+#     'U': (recip_lat_basis * np.array([0.5*gamma_bi+0.25, 1.0-gamma_bi, 0.5*gamma_bi+0.25]).T).tolist()
 # }