diff --git a/dpnegf/negf/density.py b/dpnegf/negf/density.py
index 4b06486..19e6ea2 100644
--- a/dpnegf/negf/density.py
+++ b/dpnegf/negf/density.py
@@ -267,7 +267,7 @@ def __init__(self, n_gauss=None, integrate_way:str="direct", e_grid=None):
 
     def density_integrate_Fiori(self,e_grid,kpoint,Vbias,block_tridiagonal,subblocks,integrate_way,deviceprop,
                                 device_atom_norbs,potential_at_atom,with_Dirichlet_leads,free_charge,E_ref,
-                                eta_lead=1e-5, eta_device=1e-5):
+                                eta_lead=1e-5, eta_device=1e-5, self_energy_save_path=None):
 
         for eidx, e in enumerate(self.integrate_range):
             if not with_Dirichlet_leads:
@@ -276,14 +276,13 @@ def density_integrate_Fiori(self,e_grid,kpoint,Vbias,block_tridiagonal,subblocks
                 #   to ensure the Neuamnn boundary condition in the leads.
                 deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save=False)
                 deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save=False)
+            else:
+                # For the Dirichlet leads, the self-energy of the leads is only calculated once and saved.
+                # In each iteration, load the pre-computed self-energy from disk.
+                deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save_path=self_energy_save_path)
+                deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save_path=self_energy_save_path)
+
 
-            # else:
-            #     # For the Dirichlet leads, the self-energy of the leads is only calculated once and saved.
-            #     # In each iteration, the self-energy of the leads is not updated.
-            #     deviceprop.lead_L.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save=True)
-            #     deviceprop.lead_R.self_energy(kpoint=kpoint, energy=e, eta_lead=eta_lead, save=True)
-            
-             
             deviceprop.cal_green_function(energy=e, kpoint=kpoint, block_tridiagonal=block_tridiagonal,\
                                           eta_device=eta_device,Vbias = Vbias)
 
diff --git a/dpnegf/runner/NEGF.py b/dpnegf/runner/NEGF.py
index cfcedbd..2681d29 100644
--- a/dpnegf/runner/NEGF.py
+++ b/dpnegf/runner/NEGF.py
@@ -588,12 +588,12 @@ def negf_compute(self,scf_require=False,Vbias=None):
                             self.negf_hamiltonian.subblocks = self.negf_hamiltonian.get_hs_device(only_subblocks=True)
                     
                     self.density.density_integrate_Fiori(
-                        e_grid = self.uni_grid, 
+                        e_grid = self.uni_grid,
                         kpoint=k,
                         Vbias=Vbias,
                         block_tridiagonal=self.block_tridiagonal,
                         subblocks=self.negf_hamiltonian.subblocks,
-                        integrate_way = self.density_options["integrate_way"],  
+                        integrate_way = self.density_options["integrate_way"],
                         deviceprop=self.deviceprop,
                         device_atom_norbs=self.device_atom_norbs,
                         potential_at_atom = self.potential_at_atom,
@@ -601,7 +601,8 @@ def negf_compute(self,scf_require=False,Vbias=None):
                         free_charge = self.free_charge,
                         eta_lead = self.eta_lead,
                         eta_device = self.eta_device,
-                        E_ref = self.deviceprop.E_ref
+                        E_ref = self.deviceprop.E_ref,
+                        self_energy_save_path = self.self_energy_save_path
                         )
                 else:
                     # TODO: add Ozaki support for NanoTCAD-style SCF