Create negf branch (#222)

* create poisson_grid_init.py

* add init to grid

* add find_atom_index

* add interface3D class

* add gate potential and boundary_potential

* add contents to interface3D

* rename poisson_gird_init to poisson_init

* add region definition to gate and medium

* add grid-point surface along x-axis

* add surface_grid

* add flux terms

* add func construction_poisson

* add poisson_scf.py

* add vbias term to device_preperty.py

* rename Dielectric

* add poisson_negf_scf to NEGF.py

* find potential on each atom gird through atom_index_dict

* add  update electron density for solving Poisson equation in poisson_negf_scf

* add scf break operation in poisson_negf_scf

* add negf_compute in the final part of poisson_negf_scf

* add poisson related json and read operation in NEGF

* add eps0 in unit of F/Angstrom

* delete poisson_scf.py, which has been added to NEGF.py

* clean import

* fix import issue and  str2float in gate and dieelctric definition

* add argcheck for gate and dielectric

* add grid

* update for run

* fix pyamg_solver.solve bug

* add tolerance et. to argcheck.py

* add some comments in NEGF.py and poisson_init

* update poisson_negf_scf

* update negf_compute

* update NEGF.py and poisson_init.py

* update density.py

* update device_property.py

* find the bug calculating different kpoints but get same result

* add  newK_flag and newV_flag

* add comments

* make Hartree2eV unit constant the same as that for dptb-negf

* fix self energy transpose in lead_property.py

* update NEGF.py

* update density.py

* update device_property.py

* update lead_property.py

* add notes in lead_property.py

* update negf_hamiltonian_init.py

* update a series of py files in negf, including poisson_init.py

* update recursive_green_cal.py

* update surface_green.py

* update argcheck.py

* add abs to formula powerlaw

* simplify real-space gird setting

* simplify boundary_points_get in poisson_init.py

* update interface class initialization

* change gate to dielectric

* add some notes

* seperate scipy and pyamg clearer

* transform csr to lil when attribtue values

* combine construct_Jac and construct_B to a func

* simplify the code for boundary conditions

* change sign from hLL+v*sLL to hLL-v*sLL

* remove unnecessary codes in poisson_init.py

* simplify solve_poisson_NRcycle in poisson_init.py

* add  proffiler in NEGF.py

* replace print by log.info

* get_hs_lead only when voltage changes

* move print to log

* log info update

* simplify NEGF.py

* move Fiori cal to density.py and simplify negf_compute

* add gauss integrate to Fiori method

* merge Fiori_gauss and Fiori_direct together

* test

* update density.py device_property.py lead_property.py for ksum

* update negf_hamiltonian_init.py for ksum

* update unit tests for ksum and poisson-dev

* update for ksum in poisson_dev

* add atom_cor and gird_cor consistency check

* add loginfo for k-point in NEGF.py

* add structure for tbtrans

* copy some test files from v1

* update NEGF.py init

* move AtomicData creator to negf_hamiltonian_init.py

* modify init of negf_hamiltonian_init.py

* add self.e3H in negf_hamiltonian_init.py

* add notes

* add atom_norb in hr2hk

* update NEGFHamiltonianInit initialize part for new api

* add pbc_negf in NEGFHamiltonianInit init

* modify pbc setting in device/lead

* add dptb/data/try_test.ipynb to gitignore

* add test data for negf

* fix api for device_property and density_Ozaki

* fix hamilotonian_init temporarily

* add pyamg() and scipy() in argcheck

* add soc in test_negf input

* add new model update to test_tbtrans_init.py temporarily

* comment import from tbtrans_init  temporarily

* modify test_tbtrans files for v2

* update test_negf_hamiltonian_init.py

* modify test_negf_device_property

* add h_device compare

* update negf_hamiltonian_init

* add print ldos in test_negf_device_property

* add eta_lead option in density_integrate_Fiori

* modify pbc and add dtype float 64 in negf_hamiltonian_init.py

* add TODO comment in negf_hamitonian_init

* add some comments

* rename path

* force test for negf in float64

* add test files for negf

* add task negf option in run.py

* set torch default float32 after extract H and S

* update test_negf_run.py for v2

* add sort_btd.py and split_btd.py

* add tbtrans_init in run.py

* add tbtrans_init in run.py and modify tbtrans_init.py

* add unit test for tbtrans

* add btd term in test_negf*

* fix pyamg import error

* update docstring in tbtrans_init

* rename some variables

* add BTD

* fix gitignore

* fix pyinstrument import  probelm in pytest

* fix pyinstrument import in pytest

* write code in negf_hamiltonian_init in compact form

* fix overlap Bool input probelm

* fix ldos and dos bug for BTD form

* add list for pbc setting in AtomicData_options

* fix ldos and dos  when length of up and down blocklist equals 0

* remove unnecessary n_proj_atom_pre variables

* fix energy setting problem in selfenergy calculation

* add test files for negf with overlap

* add test_negf_hamiltonian_init file

* remove unnecssary band.json

* add .gitkeep in directory

* add poisson options in argcheck

* small fix

* fix Vbias setting caused in sorting when using BTD

* update uni_test

* fix bugs and keep dptb-Poisson SCF same as that in NanoTCAD restrictly.

* update density calculation with BTD and non-BTD consistent with NanoTCAD

* update test files

* add show_blocks in BTD mode

* fix complex warning

* add DFTB+ band plot tool

* fix Fiori import error when integrate way is direct

* fix kpoints as nested_tensor form

* fix Fiori n_gauss import bug

* recover NEGF related files

* update try_bloch.ipynb

* add Bloch theorem in NEGF module to accelarate self-energy calculation

* fix bloch-fold self energy calculation bug and run successfully

* set dtype as float64 for unittest

* update try_bloch.ipynb

* support useBloch and bloch_factor in json format

* add onsite_shift in loss.py HamilLossAnalysis

* fix useBloch bug in NEGF.py

* fix sort_indices bug in get_lead_structure

* update files for bloch-self-energy in unittest

* update unittest files for hamiltonian.initialize's update

* use rmse as check standard in comparing HK_lead

* remove self energy saver temporarily

* fix pbc bug to accomdate the removing of pbc in AtomicData_options

* remove pbc output

* add overlap output in feature_to_block

* set rmse_l_HK/SK as 0 when using Block

* add block_tridiagonal in compute_density_Ozaki

* add block_tri in deviceprop.cal_green_function

* update  NEGF.py  for poisson-NEGF SCF

* add dope region in poisson-NEGF SCF

* add dope region and I-V curve calculation

* set r_max oer_max and er_max default value as float type

* change remove_bonds_nonpbc as staticmethod and add h_lead_threshold to localize the lead hamiltonian

* enlarge myConvTest in surface_green to ensure the stablization of self energy algorithm

* make correction for electro-boundary condition at Gate

* output free charge with Fiori method and seperate negf_compute from poisson_scf

* add doc to bloch

* update docstring in bloch

* add correct sign to Gate.Ef in unit of eV

.gitignore

@@ -21,6 +21,10 @@ dptb/tests/data/test_all/checkpoint/best_nnsk_b5.000_c6.615_w0.265.json
 dptb/tests/data/test_all/checkpoint/best_nnsk_b4.000_c4.000_w0.300.json
 dptb/tests/data/test_all/fancy_ones/checkpoint/best_nnsk_b4.000_c4.000_w0.300.json
 dptb/tests/data/test_negf/test_negf_run/out_negf/run_config.json
+dptb/data/try_test.ipynb
+dptb/negf/check.ipynb
+dptb/tests/data/test_negf/show.ipynb
+dptb/tests/data/test_tbtrans/show.ipynb
 run_config.json
 dptb/nnet/__pycache__/
 dptb/sktb/__pycache__/
@@ -172,3 +176,6 @@ dmypy.json
 _date.py
 _version.py
 /.idea/
+
+
+

dptb/entrypoints/run.py

@@ -9,10 +9,15 @@
 from dptb.utils.argcheck import normalize_run
 from dptb.utils.tools import j_loader
 from dptb.utils.tools import j_must_have
+from dptb.postprocess.NEGF import NEGF
+from dptb.postprocess.tbtrans_init import TBTransInputSet,sisl_installed
+
+# from dptb.postprocess.write_ham import write_ham
 from dptb.postprocess.write_block import write_block
 import torch
 import h5py
 
+
 log = logging.getLogger(__name__)
 
 def run(
@@ -87,6 +92,48 @@ def run(
                         emax=jdata["task_options"].get("emax", None))
         log.info(msg='band calculation successfully completed.')
 
+    elif task=='negf':
+
+        # try:
+        #     from pyinstrument import Profiler
+        # except ImportWarning:
+        #     log.warning(msg="pyinstrument is not installed, no profiling will be done.")
+        #     Profiler = None
+        # if Profiler is not None:
+        #     profiler = Profiler()
+        #     profiler.start()
+
+        negf = NEGF(
+            model=model,
+            AtomicData_options=jdata['AtomicData_options'],
+            structure=structure,
+            results_path=results_path,  
+            **task_options
+            )
+
+        negf.compute()
+        log.info(msg='negf calculation successfully completed.')
+
+        # if Profiler is not None:
+        #     profiler.stop()
+        #     with open(results_path+'/profile_report.html', 'w') as report_file:
+        #         report_file.write(profiler.output_html())
+
+    elif task == 'tbtrans_negf':
+        if not(sisl_installed):
+            log.error(msg="sisl is required to perform tbtrans calculation !")
+            raise RuntimeError
+        basis_dict = json.load(open(init_model))['common_options']['basis']
+        tbtrans_init = TBTransInputSet(
+            model=model, 
+            AtomicData_options=jdata['AtomicData_options'],
+            structure=structure,
+            basis_dict=basis_dict,
+            results_path=results_path,
+            **task_options)
+        tbtrans_init.hamil_get_write(write_nc=True)
+        log.info(msg='TBtrans input files are successfully generated.')
+
     elif task=='write_block':
         task = torch.load(init_model, map_location="cpu")["task"]
         block = write_block(data=struct_file, AtomicData_options=jdata['AtomicData_options'], model=model, device=jdata["device"])
@@ -95,4 +142,4 @@ def run(
             default_group = fid.create_group("0")
             for key_str, value in block.items():
                 default_group[key_str] = value.detach().cpu().numpy()
-        log.info(msg='write block successfully completed.')
+        log.info(msg='write block successfully completed.')

dptb/negf/bloch.py

@@ -0,0 +1,57 @@
+import numpy as np
+
+
+# The Bloch class in Python defines a method to unfold k points based on a given Bloch factor.
+class Bloch(object):
+
+    def __init__(self,bloch_factor) -> None:
+        '''This Python function initializes an object with a Bloch factor represented as a 3D vector.
+        
+        Parameters
+        ----------
+        bloch_factor
+            It is expected to be a list or numpy array representing a 3D vector to expand the provided
+         k points.
+        
+        '''
+
+        if  isinstance(bloch_factor,list):
+            bloch_factor = np.array(bloch_factor)
+
+        assert bloch_factor.shape[0] == 3, "kpoint should be a 3D vector"
+        self.bloch_factor = bloch_factor
+
+
+    def unfold_points(self,k:list) -> np.ndarray:
+        '''The `unfold_points` function generates expansion k points based on Bloch theorem and reshapes the
+        output into a specific format.
+        
+        Parameters
+        ----------
+        k
+            The `k` parameter in the `unfold_points` method represents the original k-point in the Brillouin zone.
+        
+        Returns
+        -------
+            The `unfold_points` method returns a reshaped array of expansion points calculated based on the
+        input parameter `k` and the bloch factor `B`.
+        
+        '''
+
+        # check k is a 3D vector
+        if isinstance(k,list):
+            assert len(k) == 3, "kpoint should be a 3D vector"
+        elif isinstance(k,np.ndarray):
+            assert k.shape[0] == 3, "kpoint should be a 3D vector"
+        else:
+            raise ValueError("k should be a list or numpy array")
+
+        # Create expansion points
+        B = self.bloch_factor
+        unfold = np.empty([B[2], B[1], B[0], 3])
+        # Use B-casting rules (much simpler)
+        unfold[:, :, :, 0] = (np.arange(B[0]).reshape(1, 1, -1) + k[0]) / B[0]
+        unfold[:, :, :, 1] = (np.arange(B[1]).reshape(1, -1, 1) + k[1]) / B[1]
+        unfold[:, :, :, 2] = (np.arange(B[2]).reshape(-1, 1, 1) + k[2]) / B[2]
+        # Back-transform shape
+        return unfold.reshape(-1, 3)

dptb/negf/density.py

@@ -140,7 +140,7 @@ def __init__(self, R, M_cut, n_gauss):
         self.R = R
         self.n_gauss = n_gauss
 
-    def integrate(self, deviceprop, kpoint, eta_lead=1e-5, eta_device=0.):
+    def integrate(self, deviceprop, kpoint, eta_lead=1e-5, eta_device=0.,Vbias=None,block_tridiagonal=False):
         '''calculates the equilibrium and non-equilibrium density matrices for a given k-point.
         
         Parameters
@@ -166,13 +166,15 @@ def integrate(self, deviceprop, kpoint, eta_lead=1e-5, eta_device=0.):
         poles = 1j* self.poles * kBT + deviceprop.lead_L.mu - deviceprop.mu # left lead expression for rho_eq
         deviceprop.lead_L.self_energy(kpoint=kpoint, energy=1j*self.R-deviceprop.mu)
         deviceprop.lead_R.self_energy(kpoint=kpoint, energy=1j*self.R-deviceprop.mu)
-        deviceprop.cal_green_function(energy=1j*self.R-deviceprop.mu, kpoint=kpoint, block_tridiagonal=False)
+        deviceprop.cal_green_function(energy=1j*self.R-deviceprop.mu, kpoint=kpoint, block_tridiagonal=block_tridiagonal,
+                                      Vbias = Vbias)
         g0 = deviceprop.grd[0]
         DM_eq = 1.0j * self.R * g0
         for i, e in enumerate(poles):
             deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
             deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
-            deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=False, eta_device=eta_device)
+            deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=block_tridiagonal, eta_device=eta_device,\
+                                          Vbias = Vbias)
             term = ((-4 * 1j * kBT) * deviceprop.grd[0] * self.residues[i]).imag
             DM_eq -= term
 
@@ -187,7 +189,7 @@ def integrate(self, deviceprop, kpoint, eta_lead=1e-5, eta_device=0.):
             for i, e in enumerate(xs):
                 deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
                 deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
-                deviceprop.cal_green_function(e=e, kpoint=kpoint, block_tridiagonal=False, eta_device=eta_device)
+                deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=block_tridiagonal, eta_device=eta_device)
                 gr_gamma_ga = torch.mm(torch.mm(deviceprop.grd[0], deviceprop.lead_R.gamma), deviceprop.grd[0].conj().T).real
                 gr_gamma_ga = gr_gamma_ga * (deviceprop.lead_R.fermi_dirac(e+deviceprop.mu) - deviceprop.lead_L.fermi_dirac(e+deviceprop.mu))
                 DM_neq = DM_neq + wlg[i] * gr_gamma_ga
@@ -225,4 +227,202 @@ def get_density_onsite(self, deviceprop, DM):
 
         onsite_density = torch.cat([torch.from_numpy(deviceprop.structure.positions), onsite_density.unsqueeze(-1)], dim=-1)
 
-        return onsite_density
+        return onsite_density
+
+
+
+class Fiori(Density):
+
+    def __init__(self, n_gauss=None):
+        super(Fiori, self).__init__()
+        self.n_gauss = n_gauss
+        self.xs = None
+        self.wlg = None
+        self.e_grid_Fiori = None
+
+    def density_integrate_Fiori(self,e_grid,kpoint,Vbias,block_tridiagonal,subblocks,integrate_way,deviceprop,
+                                device_atom_norbs,potential_at_atom,free_charge, 
+                                eta_lead=1e-5, eta_device=1e-5):
+        if integrate_way == "gauss":
+            assert self.n_gauss is not None, "n_gauss must be set in the Fiori class"
+            if self.xs is None:
+                self.xs, self.wlg = gauss_xw(xl=torch.scalar_tensor(e_grid[0]), xu=torch.scalar_tensor(e_grid[-1]), n=self.n_gauss)
+                # self.xs = self.xs.numpy();self.wlg = self.wlg.numpy()
+                self.e_grid_Fiori = e_grid
+            elif self.e_grid_Fiori[0] != e_grid[0] or self.e_grid_Fiori[-1] != e_grid[-1]:
+                self.xs, self.wlg = gauss_xw(xl=torch.scalar_tensor(e_grid[0]), xu=torch.scalar_tensor(e_grid[-1]), n=self.n_gauss)
+                # self.xs = self.xs.numpy();self.wlg = self.wlg.numpy()
+                self.e_grid_Fiori = e_grid
+            integrate_range = self.xs
+            pre_factor = self.wlg
+        elif integrate_way == "direct":
+            dE = e_grid[1] - e_grid[0]
+            integrate_range = e_grid
+            pre_factor = dE * torch.ones(len(e_grid))
+        else:
+            raise ValueError("integrate_way only supports gauss and direct in this version")
+
+
+
+        for eidx, e in enumerate(integrate_range):
+            deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
+            deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
+            deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=block_tridiagonal, eta_device=eta_device,Vbias = Vbias)
+
+            tx, ty = deviceprop.g_trans.shape
+            lx, ly = deviceprop.lead_L.se.shape
+            rx, ry = deviceprop.lead_R.se.shape
+            x0 = min(lx, tx)
+            x1 = min(rx, ty)
+
+            gammaL = torch.zeros(size=(tx, tx), dtype=torch.complex128)
+            gammaL[:x0, :x0] += deviceprop.lead_L.gamma[:x0, :x0]
+            gammaR = torch.zeros(size=(ty, ty), dtype=torch.complex128)
+            gammaR[-x1:, -x1:] += deviceprop.lead_R.gamma[-x1:, -x1:]
+
+            if not block_tridiagonal:
+                A_Rd = [torch.mm(torch.mm(deviceprop.grd[i],gammaR),deviceprop.grd[i].conj().T) for i in range(len(deviceprop.grd))]
+            else:
+                A_Rd = [torch.mm(torch.mm(deviceprop.gr_lc[i],gammaR[-x1:, -x1:]),deviceprop.gr_lc[i].conj().T) for i in range(len(deviceprop.gr_lc))]
+
+            A_Ld = [1j*(deviceprop.grd[i]-deviceprop.grd[i].conj().T)-A_Rd[i] for i in range(len(A_Rd))]
+            # the chemical potential in fermi_dirac is always set as lead_L.mu, representing the source and drain fermi level
+            gnd = [A_Ld[i]*deviceprop.lead_L.fermi_dirac(e+deviceprop.lead_L.efermi) \
+                    +A_Rd[i]*deviceprop.lead_R.fermi_dirac(e+deviceprop.lead_R.efermi) for i in range(len(A_Ld))]
+            gpd = [A_Ld[i] + A_Rd[i] - gnd[i] for i in range(len(A_Ld))]
+
+
+            # Vbias = -1 * potential_at_orb
+            for atom_index, Ei_at_atom in enumerate(-1*potential_at_atom):
+                pre_orbs = sum(device_atom_norbs[:atom_index])
+                last_orbs = pre_orbs + device_atom_norbs[atom_index]
+                # electron density
+                if e >= Ei_at_atom: 
+                    if not block_tridiagonal:
+                        free_charge[str(kpoint)][atom_index] +=\
+                            pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(gnd[0][pre_orbs:last_orbs,pre_orbs:last_orbs])
+                        # free_charge[str(kpoint)][atom_index] +=\
+                        #     pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(deviceprop.gnd[0][pre_orbs:last_orbs,pre_orbs:last_orbs])                            
+                    else:
+                        block_indexs,orb_start,orb_end = self.get_subblock_index(subblocks,atom_index,device_atom_norbs)
+                        if len(block_indexs) == 1:
+                            free_charge[str(kpoint)][atom_index] += \
+                            pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(gnd[block_indexs[0]][orb_start:orb_end,orb_start:orb_end])
+                        else:
+                            for bindex in block_indexs:
+                                if bindex == block_indexs[0]:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(gnd[bindex][orb_start:,orb_start:])
+                                elif bindex == block_indexs[-1]:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(gnd[bindex][:orb_end,:orb_end])
+                                else:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*(-1)/2/torch.pi*torch.trace(gnd[bindex])
+                # hole density
+                else:
+                    if not block_tridiagonal:                      
+                        free_charge[str(kpoint)][atom_index] +=\
+                        pre_factor[eidx]*2/2/torch.pi*torch.trace(gpd[0][pre_orbs:last_orbs,pre_orbs:last_orbs])
+                        # free_charge[str(kpoint)][atom_index] += pre_factor[eidx]*2*1/2/torch.pi*torch.trace(gpd[0][pre_orbs:last_orbs,pre_orbs:last_orbs])
+
+                    else:
+                        block_indexs,orb_start,orb_end = self.get_subblock_index(subblocks,atom_index,device_atom_norbs)
+                        if len(block_indexs) == 1:
+                            free_charge[str(kpoint)][atom_index] += \
+                            pre_factor[eidx]*2*1/2/torch.pi*torch.trace(gpd[block_indexs[0]][orb_start:orb_end,orb_start:orb_end])
+                        else:
+                            for bindex in block_indexs:
+                                if bindex == block_indexs[0]:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*1/2/torch.pi*torch.trace(gpd[bindex][orb_start:,orb_start:])
+                                elif bindex == block_indexs[-1]:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*1/2/torch.pi*torch.trace(gpd[bindex][:orb_end,:orb_end])
+                                else:
+                                    free_charge[str(kpoint)][atom_index] += \
+                                    pre_factor[eidx]*2*1/2/torch.pi*torch.trace(gpd[bindex])
+
+    def get_subblock_index(self,subblocks,atom_index,device_atom_norbs):
+        # print('atom_index:',atom_index)
+        # print('subblocks:',subblocks)
+        subblocks_cumsum = [0]+list(np.cumsum(subblocks))
+        # print('subblocks_cumsum:',subblocks_cumsum)
+        pre_orbs = sum(device_atom_norbs[:atom_index])
+        last_orbs = pre_orbs + device_atom_norbs[atom_index]
+
+        # print('pre_orbs:',pre_orbs)
+        # print('last_orbs:',last_orbs)
+
+        block_index = []
+        for i in range(len(subblocks_cumsum)-1):
+            if pre_orbs >= subblocks_cumsum[i] and last_orbs <= subblocks_cumsum[i+1]:
+                block_index.append(i)
+                orb_start = pre_orbs - subblocks_cumsum[i]
+                orb_end = last_orbs - subblocks_cumsum[i]
+                # print('1')
+                break
+            elif pre_orbs >= subblocks_cumsum[i] and pre_orbs < subblocks_cumsum[i+1] and last_orbs > subblocks_cumsum[i+1]:
+                block_index.append(i)
+                orb_start = pre_orbs - subblocks_cumsum[i]
+                for j in range(i+1,len(subblocks_cumsum)-1):
+                    block_index.append(j)
+                    if last_orbs <= subblocks_cumsum[j+1]:
+                        orb_end = last_orbs - subblocks_cumsum[j]
+                        # print('2')
+                        break
+        # print('block_index',block_index)
+        # print('orb_start',orb_start)
+        # print('orb_end',orb_end)
+        return block_index,orb_start,orb_end                 
+
+
+
+    # def density_integrate_Fiori_gauss(self,e_grid,kpoint,Vbias,deviceprop,device_atom_norbs,potential_at_atom,free_charge, eta_lead=1e-5, eta_device=1e-5):
+
+    #     if self.xs is None:
+    #         self.xs, self.wlg = gauss_xw(xl=torch.scalar_tensor(e_grid[0]), xu=torch.scalar_tensor(e_grid[-1]), n=self.n_gauss)
+    #         # self.xs = self.xs.numpy();self.wlg = self.wlg.numpy()
+    #         self.e_grid_Fiori = e_grid
+    #     elif self.e_grid_Fiori[0] != e_grid[0] or self.e_grid_Fiori[-1] != e_grid[-1]:
+    #         self.xs, self.wlg = gauss_xw(xl=torch.scalar_tensor(e_grid[0]), xu=torch.scalar_tensor(e_grid[-1]), n=self.n_gauss)
+    #         # self.xs = self.xs.numpy();self.wlg = self.wlg.numpy()
+    #         self.e_grid_Fiori = e_grid
+
+    #     for eidx, e in enumerate(self.xs):
+
+    #         deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
+    #         deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead)
+    #         deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=False, eta_device=eta_device,Vbias = Vbias)
+
+    #         tx, ty = deviceprop.g_trans.shape
+    #         lx, ly = deviceprop.lead_L.se.shape
+    #         rx, ry = deviceprop.lead_R.se.shape
+    #         x0 = min(lx, tx)
+    #         x1 = min(rx, ty)
+
+    #         gammaL = torch.zeros(size=(tx, tx), dtype=torch.complex128)
+    #         gammaL[:x0, :x0] += deviceprop.lead_L.gamma[:x0, :x0]
+    #         gammaR = torch.zeros(size=(ty, ty), dtype=torch.complex128)
+    #         gammaR[-x1:, -x1:] += deviceprop.lead_R.gamma[-x1:, -x1:]
+
+    #         A_L = torch.mm(torch.mm(deviceprop.g_trans,gammaL),deviceprop.g_trans.conj().T)
+    #         A_R = torch.mm(torch.mm(deviceprop.g_trans,gammaR),deviceprop.g_trans.conj().T)
+
+    #         # Vbias = -1 * potential_at_orb
+    #         for atom_index, Ei_at_atom in enumerate(-1*potential_at_atom):
+    #             pre_orbs = sum(device_atom_norbs[:atom_index])
+
+    #             # electron density
+    #             if e >= Ei_at_atom: 
+    #                 for j in range(device_atom_norbs[atom_index]):
+    #                     free_charge[str(kpoint)][atom_index] +=\
+    #                     self.wlg[eidx]*2*(-1)/2/torch.pi*(A_L[pre_orbs+j,pre_orbs+j]*deviceprop.lead_L.fermi_dirac(e+deviceprop.lead_L.efermi) \
+    #                                 +A_R[pre_orbs+j,pre_orbs+j]*deviceprop.lead_R.fermi_dirac(e+deviceprop.lead_R.efermi))
+
+    #             # hole density
+    #             else:
+    #                 for j in range(device_atom_norbs[atom_index]):
+    #                     free_charge[str(kpoint)][atom_index] +=\
+    #                     self.wlg[eidx]*2*1/2/torch.pi*(A_L[pre_orbs+j,pre_orbs+j]*(1-deviceprop.lead_L.fermi_dirac(e+deviceprop.lead_L.efermi)) \
+    #                                 +A_R[pre_orbs+j,pre_orbs+j]*(1-deviceprop.lead_R.fermi_dirac(e+deviceprop.lead_R.efermi)))