Merge pull request #2 from mailhexu/master

Miscellaneous

TB2J/Oiju.py

@@ -45,10 +45,6 @@ def gen_exchange_Oiju(
     kmesh=[5, 5, 5],
     emin=-12.0,
     emax=0.0,
-    height=0.2,
-    nz1=50,
-    nz2=200,
-    nz3=50,
     exclude_orbs=[],
     Rmesh=[1, 1, 1],
     description="",
@@ -93,10 +89,6 @@ def gen_exchange_Oiju(
             kmesh=kmesh,
             emin=emin,
             emax=emax,
-            height=height,
-            nz1=nz1,
-            nz2=nz2,
-            nz3=nz3,
             exclude_orbs=exclude_orbs,
             Rmesh=Rmesh,
             description=description,
@@ -123,10 +115,6 @@ def gen_exchange_Oiju(
         kmesh=[4, 4, 4],
         emin=-12.0,
         emax=0.0,
-        height=0.1,
-        nz1=50,
-        nz2=200,
-        nz3=50,
         exclude_orbs=[],
         Rmesh=[1, 1, 1],
         description="",

TB2J/Oiju_epc.py

@@ -19,10 +19,6 @@ def __init__(
         kmesh=[4, 4, 4],
         emin=-15,  # integration lower bound, relative to fermi energy
         emax=0.05,  # integration upper bound. Should be 0 (fermi energy). But DFT codes define Fermi energy in various ways.
-        height=0.5,  # the delta in the (i delta) in green's function to prevent divergence
-        nz1=150,  # grid from emin to emin+(i delta)
-        nz2=300,  # grid from emin +(i delta) to emax+(i delta)
-        nz3=150,  # grid from emax + (i delta) to emax
         exclude_orbs=[],  #
         Rmesh=[0, 0, 0],  # Rmesh.
         description="",
@@ -40,10 +36,6 @@ def __init__(
             emin=-15,  # integration lower bound, relative to fermi energy
             emax=0.05,  # integration upper bound. Should be 0 (fermi energy).
             # But DFT codes define Fermi energy in various ways.
-            height=0.5,  # the delta in the (i delta) in green's function.
-            nz1=150,  # grid from emin to emin+(i delta)
-            nz2=300,  # grid from emin +(i delta) to emax+(i delta)
-            nz3=150,  # grid from emax + (i delta) to emax
             exclude_orbs=[],  #
             Rmesh=[0, 0, 0],  # Rmesh.
             description="",

TB2J/__init__.py

@@ -1 +1 @@
-__version__ = "0.8.2.4"
+__version__ = "0.8.2.8"

TB2J/abacus/abacus_wrapper.py

@@ -18,8 +18,7 @@ class AbacusWrapper(AbstractTB):
     def __init__(self, HR, SR, Rlist, nbasis, nspin=1):
         self.R2kfactor = -2j * np.pi
         self.is_orthogonal = False
-        self.is_siesta = False
-        self._name = "Abacus"
+        self._name = "ABACUS"
         self.HR = HR
         self.SR = SR
         self.Rlist = Rlist

TB2J/exchange.py

@@ -26,10 +26,6 @@ def __init__(
         emax=0.05,
         nz=100,
         # the delta in the (i delta) in green's function to prevent divergence
-        height=0.5,
-        nz1=150,  # grid from emin to emin+(i delta)
-        nz2=300,  # grid from emin +(i delta) to emax+(i delta)
-        nz3=150,  # grid from emax + (i delta) to emax
         exclude_orbs=[],  #
         ne=None,  # number of electrons in Wannier function.
         Rcut=None,  # Rcut.
@@ -44,12 +40,6 @@ def __init__(
         self.emin = emin
         self.emax = emax
         self.nz = nz
-        self.height = height
-        self.nz1 = nz1
-        self.nz2 = nz2
-        self.nz3 = nz3
-        if nz is None:
-            self.nz = nz1 + nz2 + nz3
         self.Rcut = Rcut
         self.basis = basis
         self.magnetic_elements = magnetic_elements
@@ -81,10 +71,6 @@ def __init__(
         emax=0.05,
         nz=100,
         # the delta in the (i delta) in green's function to prevent divergence
-        height=0.5,
-        nz1=150,  # grid from emin to emin+(i delta)
-        nz2=300,  # grid from emin +(i delta) to emax+(i delta)
-        nz3=150,  # grid from emax + (i delta) to emax
         exclude_orbs=[],  #
         ne=None,  # number of electrons in Wannier function.
         Rcut=None,  # Rcut.
@@ -105,10 +91,6 @@ def __init__(
             emin=emin,
             emax=emax,
             nz=nz,
-            height=height,
-            nz1=nz1,
-            nz2=nz2,
-            nz3=nz3,
             exclude_orbs=exclude_orbs,
             ne=ne,
             Rcut=Rcut,
@@ -144,6 +126,7 @@ def _prepare_Jorb_file(self):
 
     def _adjust_emin(self):
         self.emin = self.G.find_energy_ingap(rbound=self.efermi - 5.0) - self.efermi
+        print(f"A gap is found at {self.emin}, set emin to it.")
 
     def set_tbmodels(self, tbmodels):
         pass
@@ -163,11 +146,11 @@ def _prepare_elist(self, method="legendre"):
          emin --1-> emin + 1j*height --2-> emax+1j*height --3-> emax
         """
         self.contour = Contour(self.emin, self.emax)
-        if method.lower() == "rectangle":
-            self.contour.build_path_rectangle(
-                height=self.height, nz1=self.nz1, nz2=self.nz2, nz3=self.nz3
-            )
-        elif method.lower() == "semicircle":
+        # if method.lower() == "rectangle":
+        #    self.contour.build_path_rectangle(
+        #        height=self.height, nz1=self.nz1, nz2=self.nz2, nz3=self.nz3
+        #    )
+        if method.lower() == "semicircle":
             self.contour.build_path_semicircle(npoints=self.nz, endpoint=True)
         elif method.lower() == "legendre":
             self.contour.build_path_legendre(npoints=self.nz, endpoint=True)
@@ -257,7 +240,7 @@ def _prepare_orb_mmat(self):
         self.mmats = {}
         self.orbital_names = {}
         self.norb_reduced = {}
-        if self.backend_name == "SIESTA":
+        if self.backend_name.upper() == "SIESTA":
             syms = self.atoms.get_chemical_symbols()
             for iatom, orbs in self.labels.items():
                 if (self.include_orbs is not None) and syms[iatom] in self.include_orbs:
@@ -322,7 +305,7 @@ def simplify_orbital_contributions(self, Jorbij, iatom, jatom):
         """
         sum up the contribution of all the orbitals with same (n,l,m)
         """
-        if self.backend_name == "SIESTA":
+        if self.backend_name.upper() == "SIESTA":
             mmat_i = self.mmats[iatom]
             mmat_j = self.mmats[jatom]
             Jorbij = mmat_i.T @ Jorbij @ mmat_j
@@ -361,7 +344,8 @@ def set_tbmodels(self, tbmodels):
         )
         self.norb = self.G.norb
         self.nbasis = self.G.nbasis
-        self.rho = np.zeros((self.nbasis, self.nbasis), dtype=complex)
+        # self.rho = np.zeros((self.nbasis, self.nbasis), dtype=complex)
+        self.rho = self.G.get_density_matrix().real
         self.A_ijR_list = defaultdict(lambda: [])
         self.A_ijR = defaultdict(lambda: np.zeros((4, 4), dtype=complex))
         self.A_ijR_orb = dict()
@@ -594,22 +578,22 @@ def A_to_Jtensor(self):
                 B = np.imag(val[3, 3])
                 self.B[keyspin] = Jprime, B
 
-    def get_N_e(self, GR, de):
-        """
-        calcualte density matrix for all R,i, j
-        """
-        self.N = defaultdict(lambda: 0.0)
-        for R, G in GR.items():
-            self.N[R] += -1.0 / np.pi * np.imag(G * de)
+    # def get_N_e(self, GR, de):
+    #    """
+    #    calcualte density matrix for all R,i, j
+    #    """
+    #    self.N = defaultdict(lambda: 0.0)
+    #    for R, G in GR.items():
+    #        self.N[R] += -1.0 / np.pi * np.imag(G * de)
 
-    def get_rho_e(self, rhoR):
-        """add component to density matrix from a green's function
-        :param GR: Green's funciton in real space.
-        """
-        return -1.0 / np.pi * rhoR[0, 0, 0]
+    # def get_rho_e(self, rhoR):
+    #    """add component to density matrix from a green's function
+    #    :param GR: Green's funciton in real space.
+    #    """
+    #    return -1.0 / np.pi * rhoR[0, 0, 0]
 
-    def get_total_charges(self):
-        return np.sum(np.imag(np.diag(self.rho)))
+    # def get_total_charges(self):
+    #    return np.sum(np.imag(np.diag(self.rho)))
 
     def get_rho_atom(self):
         """
@@ -622,9 +606,9 @@ def get_rho_atom(self):
             iorb = self.iorb(iatom)
             tmp = self.rho[np.ix_(iorb, iorb)]
             # *2 because there is a 1/2 in the paui_block_all function
-            rho[iatom] = np.array([np.trace(x) * 2 for x in pauli_block_all(tmp)])
-            self.charges[iatom] = np.imag(rho[iatom][0])
-            self.spinat[iatom, :] = np.imag(rho[iatom][1:])
+            rho[iatom] = np.array([np.trace(x) * 2 for x in pauli_block_all(tmp)]).real
+            self.charges[iatom] = rho[iatom][0]
+            self.spinat[iatom, :] = rho[iatom][1:]
         self.rho_dict = rho
         return self.rho_dict
 
@@ -675,7 +659,7 @@ def calculate_DMI_NJT(self):
         self.Ddict_NJT = Ddict_NJT
         return Ddict_NJT
 
-    def integrate(self, rhoRs, AijRs, AijRs_orb=None, method="simpson"):
+    def integrate(self, AijRs, AijRs_orb=None, method="simpson"):
         """
         AijRs: a list of AijR,
         wherer AijR: array of ((nR, n, n, 4,4), dtype=complex)
@@ -685,7 +669,7 @@ def integrate(self, rhoRs, AijRs, AijRs_orb=None, method="simpson"):
         elif method == "simpson":
             integrate = simpson_nonuniform
 
-        self.rho = integrate(self.contour.path, rhoRs)
+        # self.rho = integrate(self.contour.path, rhoRs)
         for iR, R in enumerate(self.R_ijatom_dict):
             for iatom, jatom in self.R_ijatom_dict[R]:
                 f = AijRs[(R, iatom, jatom)]
@@ -695,10 +679,10 @@ def integrate(self, rhoRs, AijRs, AijRs_orb=None, method="simpson"):
                         self.contour.path, AijRs_orb[(R, iatom, jatom)]
                     )
 
-    def get_AijR_rhoR(self, e):
-        GR, rhoR = self.G.get_GR(self.short_Rlist, energy=e, get_rho=True)
+    def get_AijR(self, e):
+        GR = self.G.get_GR(self.short_Rlist, energy=e, get_rho=False)
         AijR, AijR_orb = self.get_all_A(GR)
-        return AijR, AijR_orb, self.get_rho_e(rhoR)
+        return AijR, AijR_orb
 
     def save_AijR(self, AijRs, fname):
         result = dict(path=self.contour.path, AijRs=AijRs)
@@ -716,7 +700,6 @@ def calculate_all(self):
         """
         print("Green's function Calculation started.")
 
-        rhoRs = []
         AijRs = {}
 
         AijRs_orb = {}
@@ -725,11 +708,9 @@ def calculate_all(self):
 
         npole = len(self.contour.path)
         if self.np > 1:
-            # executor = ProcessPool(nodes=self.np)
-            # results = executor.map(self.get_AijR_rhoR, self.contour.path)
-            results = p_map(self.get_AijR_rhoR, self.contour.path, num_cpus=self.np)
+            results = p_map(self.get_AijR, self.contour.path, num_cpus=self.np)
         else:
-            results = map(self.get_AijR_rhoR, tqdm(self.contour.path, total=npole))
+            results = map(self.get_AijR, tqdm(self.contour.path, total=npole))
 
         for i, result in enumerate(results):
             for iR, R in enumerate(self.R_ijatom_dict):
@@ -749,16 +730,14 @@ def calculate_all(self):
                             AijRs_orb[(R, iatom, jatom)].append(
                                 result[1][R, iatom, jatom]
                             )
-            rhoRs.append(result[2])
         if self.np > 1:
             # executor.close()
             # executor.join()
             # executor.clear()
             pass
 
         # self.save_AijRs(AijRs)
-        self.integrate(rhoRs, AijRs, AijRs_orb)
-
+        self.integrate(AijRs, AijRs_orb)
         self.get_rho_atom()
         self.A_to_Jtensor()
         self.A_to_Jtensor_orb()

TB2J/exchangeCL2.py

@@ -42,10 +42,10 @@ def set_tbmodels(self, tbmodels):
             )
         self.norb = self.Gup.norb
         self.nbasis = self.Gup.nbasis + self.Gdn.nbasis
-        self.rho_up_list = []
-        self.rho_dn_list = []
-        self.rho_up = np.zeros((self.norb, self.norb), dtype=float)
-        self.rho_dn = np.zeros((self.norb, self.norb), dtype=float)
+        # self.rho_up_list = []
+        # self.rho_dn_list = []
+        self.rho_up = self.Gup.get_density_matrix()
+        self.rho_dn = self.Gdn.get_density_matrix()
         self.Jorb_list = defaultdict(lambda: [])
         self.JJ_list = defaultdict(lambda: [])
         self.JJ = defaultdict(lambda: 0.0j)
@@ -225,8 +225,6 @@ def integrate(self, method="simpson"):
             integrate = trapezoidal_nonuniform
         elif method == "simpson":
             integrate = simpson_nonuniform
-        self.rho_up = np.imag(integrate(self.contour.path, self.rho_up_list))
-        self.rho_dn = np.imag(integrate(self.contour.path, self.rho_dn_list))
         for R, ijpairs in self.R_ijatom_dict.items():
             for iatom, jatom in ijpairs:
                 self.Jorb[(R, iatom, jatom)] = integrate(
@@ -236,12 +234,11 @@ def integrate(self, method="simpson"):
                     self.contour.path, self.JJ_list[(R, iatom, jatom)]
                 )
 
-    def get_AijR_rhoR(self, e):
-        GR_up, rho_up = self.Gup.get_GR(self.short_Rlist, energy=e, get_rho=True)
-        GR_dn, rho_dn = self.Gdn.get_GR(self.short_Rlist, energy=e, get_rho=True)
-        rup, rdn = self.get_rho_e(rho_up, rho_dn)
+    def get_AijR(self, e):
+        GR_up = self.Gup.get_GR(self.short_Rlist, energy=e, get_rho=False)
+        GR_dn = self.Gdn.get_GR(self.short_Rlist, energy=e, get_rho=False)
         Jorb_list, JJ_list = self.get_all_A(GR_up, GR_dn)
-        return rup, rdn, Jorb_list, JJ_list
+        return Jorb_list, JJ_list
 
     def calculate_all(self):
         """
@@ -251,15 +248,13 @@ def calculate_all(self):
 
         npole = len(self.contour.path)
         if self.np == 1:
-            results = map(self.get_AijR_rhoR, tqdm(self.contour.path, total=npole))
+            results = map(self.get_AijR, tqdm(self.contour.path, total=npole))
         else:
             # pool = ProcessPool(nodes=self.np)
             # results = pool.map(self.get_AijR_rhoR ,self.contour.path)
-            results = p_map(self.get_AijR_rhoR, self.contour.path, num_cpus=self.np)
+            results = p_map(self.get_AijR, self.contour.path, num_cpus=self.np)
         for i, result in enumerate(results):
-            rup, rdn, Jorb_list, JJ_list = result
-            self.rho_up_list.append(rup)
-            self.rho_dn_list.append(rdn)
+            Jorb_list, JJ_list = result
             for iR, R in enumerate(self.R_ijatom_dict):
                 for iatom, jatom in self.R_ijatom_dict[R]:
                     key = (R, iatom, jatom)

TB2J/exchange_pert.py

@@ -194,7 +194,7 @@ def calculate_all(self):
             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)
-            self.get_rho_e(GR, de)
+            # self.get_rho_e(GR, de)
             self.get_all_A(GR, dGdx, de)
             if self.calc_NJt:
                 self.get_N_e(GR, de)

TB2J/green.py

@@ -225,22 +225,29 @@ def get_Sk(self, ik):
 
     def get_density_matrix(self):
         rho = np.zeros((self.nbasis, self.nbasis), dtype=complex)
-        for ik, _ in enumerate(self.kpts):
-            rho += (
-                (self.get_evecs(ik) * fermi(self.evals[ik], self.efermi))
-                @ self.get_evecs(ik).T.conj()
-                * self.kweights[ik]
-            )
+        if self.is_orthogonal:
+            for ik, _ in enumerate(self.kpts):
+                rho += (
+                    (self.get_evecs(ik) * fermi(self.evals[ik], self.efermi))
+                    @ self.get_evecs(ik).T.conj()
+                    * self.kweights[ik]
+                )
+        else:
+            for ik, _ in enumerate(self.kpts):
+                rho += (
+                    (self.get_evecs(ik) * fermi(self.evals[ik], self.efermi))
+                    @ self.get_evecs(ik).T.conj()
+                    @ self.get_Sk(ik)
+                    * self.kweights[ik]
+                )
+
         return rho
 
     def get_rho_R(self, Rlist):
         nR = len(Rlist)
         rho_R = np.zeros((nR, self.nbasis, self.nbasis), dtype=complex)
         for ik, kpt in enumerate(self.kpts):
             evec = self.get_evecs(ik)
-            # rhok=(evec * fermi(self.evals[ik], self.efermi)
-            #        ) @ evec.T.conj()
-            # print(fermi(self.evals[ik] , self.efermi))
             rhok = np.einsum(
                 "ib,b, bj-> ij", evec, fermi(self.evals[ik], self.efermi), evec.conj().T
             )