diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml
index 7f453c0..bd3d9d0 100644
--- a/.github/workflows/python-app.yml
+++ b/.github/workflows/python-app.yml
@@ -14,9 +14,7 @@ permissions:
 
 jobs:
   build:
-
     runs-on: ubuntu-latest
-
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python 3.10
diff --git a/TB2J/Oiju_epc.py b/TB2J/Oiju_epc.py
index 7426593..ff850e3 100644
--- a/TB2J/Oiju_epc.py
+++ b/TB2J/Oiju_epc.py
@@ -57,6 +57,8 @@ def __init__(
         # prepare dDelta
         self.calc_dDelta()
 
+        self._kmap = None
+
     def prepare_epc_wann(self):
         """
         prepare EPC in electron wannier function representation
@@ -64,6 +66,12 @@ def prepare_epc_wann(self):
         self.EPCmat_wann_up = self.EPCmat_up.to_wann(self.Umat_up)
         self.EPCmat_wann_dn = self.EPCmat_dn.to_wann(self.Umat_dn)
 
+    def get_iq(self, q):
+        if self._kmap is None:
+            self._kmap = {}
+        for i, q in enumerate(self.klist):
+            self._kmap[tuple(q)] = i
+
     def calc_dDelta(self):
         """
         calculate $\delta \Delta$.
@@ -77,7 +85,8 @@ def calc_dDelta(self):
         """
         self.dDelta = np.zeros((self.nphon, self.nwann, self.nwann))
         # iq, iv, iR
-        iq0 = iqlist[(0, 0, 0)]
+        # iq0 = self.iqlist[(0, 0, 0)]
+        iq0 = self.get_iq((0, 0, 0))
 
         for iv in range(self.nphon):
             for ik, k in enumerate(self.klist):
diff --git a/TB2J/exchangeCL2.py b/TB2J/exchangeCL2.py
index b16d712..4e96b94 100644
--- a/TB2J/exchangeCL2.py
+++ b/TB2J/exchangeCL2.py
@@ -117,17 +117,17 @@ def get_A_ijR(self, Gup, Gdn, iatom, jatom):
                     "ij, ji-> ij", np.matmul(Deltai, Gij_up), np.matmul(Deltaj, Gji_dn)
                 )
 
-                if self.biquadratic:
-                    A = np.einsum(
-                        "ij, ji-> ij",
-                        np.matmul(Deltai, Gij_up),
-                        np.matmul(Deltaj, Gji_up),
-                    )
-                    C = np.einsum(
-                        "ij, ji-> ij",
-                        np.matmul(Deltai, Gij_down),
-                        np.matmul(Deltaj, Gji_down),
-                    )
+                # if self.biquadratic:
+                #    A = np.einsum(
+                #        "ij, ji-> ij",
+                #        np.matmul(Deltai, Gij_up),
+                #        np.matmul(Deltaj, Gji_up),
+                #    )
+                #    C = np.einsum(
+                #        "ij, ji-> ij",
+                #        np.matmul(Deltai, Gij_down),
+                #        np.matmul(Deltaj, Gji_down),
+                #    )
                 tmp = np.sum(t)
                 self.Jorb_list[(R, iatom, jatom)].append(t / (4.0 * np.pi))
                 self.JJ_list[(R, iatom, jatom)].append(tmp / (4.0 * np.pi))
diff --git a/TB2J/exchange_pert.py b/TB2J/exchange_pert.py
index 15c7d05..86ea877 100644
--- a/TB2J/exchange_pert.py
+++ b/TB2J/exchange_pert.py
@@ -178,19 +178,19 @@ def calculate_all(self):
         """
         print("Green's function Calculation started.")
 
-        widgets = [
-            " [",
-            progressbar.Timer(),
-            "] ",
-            progressbar.Bar(),
-            " (",
-            progressbar.ETA(),
-            ") ",
-        ]
-        bar = progressbar.ProgressBar(maxval=self.contour.npoints, widgets=widgets)
-        bar.start()
+        # widgets = [
+        #    " [",
+        #    progressbar.Timer(),
+        #    "] ",
+        #    progressbar.Bar(),
+        #    " (",
+        #    progressbar.ETA(),
+        #    ") ",
+        # ]
+        # bar = progressbar.ProgressBar(maxval=self.contour.npoints, widgets=widgets)
+        # bar.start()
         for ie in range(self.contour.npoints):
-            bar.update(ie)
+            # bar.update(ie)
             e = self.contour.elist[ie]
             de = self.contour.de[ie]
             GR, dGdx = self.G.get_GR_and_dGRdx(self.Rlist, energy=e, dHdx=self.dHdx)
@@ -203,7 +203,7 @@ def calculate_all(self):
         self.A_to_Jtensor()
         if self.calc_NJt:
             self.calculate_DMI_NJT()
-        bar.finish()
+        # bar.finish()
 
     def write_output(self, path="TB2J_results"):
         self._prepare_index_spin()
diff --git a/TB2J/exchange_qspace.py b/TB2J/exchange_qspace.py
index c95e031..61a120b 100644
--- a/TB2J/exchange_qspace.py
+++ b/TB2J/exchange_qspace.py
@@ -90,8 +90,8 @@ def get_all_A(self):
                     # Kq= Gk_up[ik] @ self.Delta @ Gk_dn[ikq]
                     # for ik, ikq in enumerate(range(self.nkpts)):
                     Guk = Gk_up[ik, :, :]
-                    # Gdk = Gk_dn[ik, :, :]
-                    # Gukq = Gk_up[ikq, :, :]
+                    Gdk = Gk_dn[ik, :, :]
+                    Gukq = Gk_up[ikq, :, :]
                     Gdkq = Gk_dn[ikq, :, :]
                     for i, iatom in enumerate(self.ind_mag_atoms):
                         Deltai = self.get_Delta(iatom)
diff --git a/TB2J/gpaw_wrapper.py b/TB2J/gpaw_wrapper.py
index fc82e6b..d8ca24b 100644
--- a/TB2J/gpaw_wrapper.py
+++ b/TB2J/gpaw_wrapper.py
@@ -61,7 +61,7 @@ def __init__(self, atoms=None, calc=None, gpw_fname=None):
     #    return eigh(self.hamk(k), self.S(k))
 
     def get_kpts(self):
-        return calc.get_ibz_k_points()
+        return self.calc.get_ibz_k_points()
 
     def HS_and_eigen(self, kpts=None, convention=2):
         if kpts is not None:
@@ -135,7 +135,7 @@ def __init__(self, calc=None, atoms=None, gpw_fname=None, pickle_fname=None):
             self.H_NMM, self.S_NMM = tb.h_and_s()
             self.Rlist = tb.R_cN.T
         else:
-            with open(fname, "rb") as myfile:
+            with open(gpw_fname, "rb") as myfile:
                 self.H_NMM, self.S_NMM, self.Rlist = pickle.load(myfile)
         self.nR, self.nbasis, _ = self.H_NMM.shape
         self.positions = np.zeros((self.nbasis, 3))
diff --git a/TB2J/mathutils.py b/TB2J/mathutils.py
new file mode 100644
index 0000000..bee9ccc
--- /dev/null
+++ b/TB2J/mathutils.py
@@ -0,0 +1,12 @@
+import numpy as np
+from scipy.linalg import inv, eigh
+
+
+def Lowdin(S):
+    """
+    Calculate S^(-1/2).
+    Which is used in lowind's symmetric orthonormalization.
+    psi_prime = S^(-1/2) psi
+    """
+    eigval, eigvec = eigh(S)
+    return eigvec @ np.diag(np.sqrt(1.0 / eigval)) @ (eigvec.T.conj())
diff --git a/TB2J/sisl_wrapper.py b/TB2J/sisl_wrapper.py
index ad0f69b..fe695c6 100644
--- a/TB2J/sisl_wrapper.py
+++ b/TB2J/sisl_wrapper.py
@@ -5,6 +5,7 @@
 from collections import defaultdict
 from scipy.linalg import eigh
 from TB2J.myTB import AbstractTB
+from TB2J.mathutils import Lowdin
 
 
 class SislWrapper(AbstractTB):
@@ -199,6 +200,9 @@ def HS_and_eigen(self, kpts, convention=2):
                     shape=(nkpts, self.nbasis, self.nbasis),
                     dtype=complex,
                 )
+            Hk = self.Hk(k, convention=convention)
+            Sk = self.Sk(k, convention=convention)
+            evalue, evec = eigh(Hk, Sk)
             self.H[ik] = Hk
             self.S[ik] = Sk
             self.evals[ik] = evalue
diff --git a/TB2J/spinham/supercell.py b/TB2J/spinham/supercell.py
index bf87007..f214642 100644
--- a/TB2J/spinham/supercell.py
+++ b/TB2J/spinham/supercell.py
@@ -319,7 +319,7 @@ def map_to_primitive(atoms, primitive_atoms, offset=(0, 0, 0)):
 def test():
     sc_mat = np.diag([1, 1, 2])
     # sc_mat[0, 1] = 2
-    spm = supercell_maker(sc_matrix=sc_mat)
+    spm = SupercellMaker(sc_matrix=sc_mat)
     print(spm.sc_cell([1, 1, 1]))
     print(spm.sc_pos([[0.5, 1, 1]]))
     print(spm.sc_trans_invariant(["Fe"]))
diff --git a/TB2J/utils.py b/TB2J/utils.py
index 10efc56..321a2ac 100644
--- a/TB2J/utils.py
+++ b/TB2J/utils.py
@@ -1,4 +1,5 @@
 from collections import OrderedDict, defaultdict
+import copy
 import numpy as np
 from pathlib import Path