add tutorial 2 without large self-energy files

.gitignore

@@ -27,6 +27,10 @@ dpnegf/data/try_test.ipynb
 dpnegf/negf/check.ipynb
 dpnegf/tests/data/test_negf/show.ipynb
 dpnegf/tests/data/test_tbtrans/show.ipynb
+# tutorial files
+dpnegf/examples/graphene/negf_output_k20/self_energy/self_energy_leadL.h5
+dpnegf/examples/graphene/negf_output_k20/self_energy/self_energy_leadR.h5
+
 run_config.json
 dpnegf/nnet/__pycache__/
 dpnegf/sktb/__pycache__/

dpnegf/runner/NEGF.py

@@ -6,7 +6,7 @@
 from dpnegf.utils.elec_struc_cal import ElecStruCal
 from dpnegf.negf.density import Ozaki,Fiori
 from dpnegf.negf.device_property import DeviceProperty
-from dpnegf.negf.lead_property import LeadProperty, compute_all_self_energy
+from dpnegf.negf.lead_property import LeadProperty, compute_all_self_energy, _has_saved_self_energy
 from dpnegf.negf.negf_utils import is_fully_covered
 import ase
 from dpnegf.utils.constants import Boltzmann, eV2J
@@ -21,7 +21,6 @@
 from pyinstrument import Profiler
 import os
 
-
 log = logging.getLogger(__name__)
 
 
@@ -49,7 +48,7 @@ def __init__(self,
                 sgf_solver: str,
                 e_fermi: float=None,
                 use_saved_HS: bool=False, saved_HS_path: str=None,
-                self_energy_save: bool=False, self_energy_save_path: str=None, se_info_display: bool=False,
+                use_saved_se: bool=False, self_energy_save_path: str=None, se_info_display: bool=False,
                 out_tc: bool=False,out_dos: bool=False,out_density: bool=False,out_potential: bool=False,
                 out_current: bool=False,out_current_nscf: bool=False,out_ldos: bool=False,out_lcurrent: bool=False,
                 results_path: Optional[str]=None,
@@ -80,9 +79,9 @@ def __init__(self,
         self.saved_HS_path = saved_HS_path
 
         self.sgf_solver = sgf_solver
-        self.self_energy_save = self_energy_save
-        self.self_energy_save_path = self_energy_save_path
-        self.se_info_display = se_info_display
+        self.use_saved_se = use_saved_se # whether to use the saved self-energy or not
+        self.self_energy_save_path = self_energy_save_path # The directory to save the self-energy or for saved self-energy
+        self.se_info_display = se_info_display # whether to display the self-energy information after calculation
         self.pbc = self.stru_options["pbc"]
 
         if  self.stru_options["lead_L"]["useBloch"] or self.stru_options["lead_R"]["useBloch"]:
@@ -133,18 +132,8 @@ def __init__(self,
             else:
                 assert self.stru_options[lead_tag]["voltage"] == 0, f"{lead_tag} voltage should be 0 in non-scf calculation"
 
-        if AtomicData_options is None:
-            from dptb.utils.argcheck import get_cutoffs_from_model_options
-            # get the cutoffs from model options
-            r_max, er_max, oer_max  = get_cutoffs_from_model_options(model.model_options)
-            AtomicData_options = {'r_max': r_max, 'er_max': er_max, 'oer_max': oer_max}
-        else:
-            log.warning(msg="AtomicData_options is extracted from input file. " \
-                            "This may be not consistent with the model options. " \
-                            "Please be careful and check the cutoffs.")
-        formatted = json.dumps(AtomicData_options, indent=4)
-        indented = '\n'.join(' ' * 15 + line for line in formatted.splitlines())
-        log.info("The AtomicData_options is:\n%s", indented)
+        # preparing AtomicData_options, including cutoffs
+        AtomicData_options = self.update_atomicdata_options(model,AtomicData_options)
 
         # computing the hamiltonian
         self.negf_hamiltonian = NEGFHamiltonianInit(model=model,
@@ -403,13 +392,13 @@ def compute(self):
             self.negf_compute(scf_require=False,Vbias=self.potential_at_orb)
 
         else:
-            # profiler = Profiler()
-            # profiler.start() 
+            profiler = Profiler()
+            profiler.start() 
             self.negf_compute(scf_require=False,Vbias=None)
-            # profiler.stop()
-            # output_path = os.path.join(self.results_path, "profile_report.html")
-            # with open(output_path, 'w') as report_file:
-                # report_file.write(profiler.output_html())
+            profiler.stop()
+            output_path = os.path.join(self.results_path, "profile_report.html")
+            with open(output_path, 'w') as report_file:
+                report_file.write(profiler.output_html())
 
     def poisson_negf_scf(self,interface_poisson,atom_gridpoint_index,err=1e-6,max_iter=1000,
                          mix_method:str='linear', mix_rate:float=0.3, tolerance:float=1e-7,Gaussian_sigma:float=3.0):
@@ -513,35 +502,57 @@ def poisson_negf_scf(self,interface_poisson,atom_gridpoint_index,err=1e-6,max_it
         # if iter_count <= max_iter: 
         #     profiler.stop()
         #     with open('profile_report.html', 'w') as report_file:
-        #         report_file.write(profiler.output_html())
-
+        #         report_file.write(profiler.output_html())、
+
+    def prepare_self_energy(self, scf_require: bool) -> None:
+        """
+        Prepares the self-energy for the NEGF calculation.
+
+        Depending on the calculation settings, this method either loads previously saved self-energy data
+        or computes and saves new self-energy values for the device leads. The computation method varies
+        based on whether self-consistent field (SCF) calculations are required and whether Dirichlet boundary
+        conditions are applied to the leads.
+
+        Parameters:
+        ----------
+            scf_require (bool): Indicates whether SCF calculations are required.
+        """
+        # self energy calculation
+        log.info(msg="------Self-energy calculation------")
+        if  self.self_energy_save_path is None:
+            self.self_energy_save_path = os.path.join(self.results_path, "self_energy") 
+        os.makedirs(self.self_energy_save_path, exist_ok=True)
+
+        if self.use_saved_se:
+            assert _has_saved_self_energy(self.self_energy_save_path), "No saved self-energy found in {}".format(self.self_energy_save_path)
+            log.info(msg="Using saved self-energy from {}".format(self.self_energy_save_path))
+            log.info(msg="Ensure the saved self-energy is consistent with the current calculation setting!")
+        else:
+            log.info(msg="Calculating self-energy and saving to {}".format(self.self_energy_save_path))
+            if scf_require and self.poisson_options["with_Dirichlet_leads"]:
+                # 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.
+                # for ik, k in enumerate(self.kpoints):
+                #     for e in self.density.integrate_range:
+                #         self.deviceprop.lead_L.self_energy(kpoint=k, energy=e, eta_lead=self.eta_lead, save=True)
+                #         self.deviceprop.lead_R.self_energy(kpoint=k, energy=e, eta_lead=self.eta_lead, save=True)
+                compute_all_self_energy(self.eta_lead, self.deviceprop.lead_L, self.deviceprop.lead_R,
+                                        self.kpoints, self.density.integrate_range, self.self_energy_save_path)
+            elif not self.scf:
+                # In non-scf case, the self-energy of the leads is calculated for each energy point in the energy grid.
+                compute_all_self_energy(self.eta_lead, self.deviceprop.lead_L, self.deviceprop.lead_R,
+                                        self.kpoints, self.uni_grid, self.self_energy_save_path)
+        log.info(msg="-----------------------------------\n")
+
+
 
     def negf_compute(self,scf_require=False,Vbias=None):
 
         assert scf_require is not None, "scf_require should be set to True or False"
         self.out['k']=[];self.out['wk']=[]
         if hasattr(self, "uni_grid"): self.out["uni_grid"] = self.uni_grid
 
-        # self energy calculation
-        log.info(msg="------Self-energy calculation------")
-        selfen_parent_dir = os.path.join(self.results_path,"self_energy")
-        if not os.path.exists(selfen_parent_dir): 
-            os.makedirs(selfen_parent_dir)
-        if scf_require and self.poisson_options["with_Dirichlet_leads"]:
-            # 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.
-            # for ik, k in enumerate(self.kpoints):
-            #     for e in self.density.integrate_range:
-            #         self.deviceprop.lead_L.self_energy(kpoint=k, energy=e, eta_lead=self.eta_lead, save=True)
-            #         self.deviceprop.lead_R.self_energy(kpoint=k, energy=e, eta_lead=self.eta_lead, save=True)
-            compute_all_self_energy(self.eta_lead, self.deviceprop.lead_L, self.deviceprop.lead_R,
-                                    self.kpoints, self.density.integrate_range)
-        elif not self.scf:
-            # In non-scf case, the self-energy of the leads is calculated for each energy point in the energy grid.
-            compute_all_self_energy(self.eta_lead, self.deviceprop.lead_L, self.deviceprop.lead_R,
-                                    self.kpoints, self.uni_grid)
-        log.info(msg="-----------------------------------\n")
-
+        self.prepare_self_energy(scf_require)
 
         for ik, k in enumerate(self.kpoints):
 
@@ -627,7 +638,6 @@ def negf_compute(self,scf_require=False,Vbias=None):
                                         kpoint=k, 
                                         eta_lead=self.eta_lead,
                                         method=self.sgf_solver,
-                                        save=self.self_energy_save,
                                         save_path=self.self_energy_save_path,
                                         se_info_display=self.se_info_display
                                         )
@@ -641,7 +651,6 @@ def negf_compute(self,scf_require=False,Vbias=None):
                                         kpoint=k, 
                                         eta_lead=self.eta_lead,
                                         method=self.sgf_solver,
-                                        save=self.self_energy_save,
                                         save_path=self.self_energy_save_path,
                                         se_info_display=self.se_info_display
                                         )                                
@@ -736,7 +745,6 @@ def negf_compute(self,scf_require=False,Vbias=None):
                                         kpoint=k, 
                                         eta_lead=self.eta_lead,
                                         method=self.sgf_solver,
-                                        save=self.self_energy_save,
                                         save_path=self.self_energy_save_path,
                                         se_info_display=self.se_info_display
                                         )
@@ -815,6 +823,42 @@ def get_nel_atom_lead(self, struct_leads, charge:float=None):
                 nel_atom_lead[lead_tag] = {elem: nel_atom_lead[lead_tag][elem] + charge[lead_tag] for elem in nel_atom_lead[lead_tag]}
 
         return nel_atom_lead  
+
+    @staticmethod
+    def update_atomicdata_options(model,AtomicData_options: dict=None) -> dict:
+        """
+        Updates or initializes the AtomicData_options dictionary based on the provided model.
+
+        If AtomicData_options is not provided, it extracts cutoff values from the model's options
+        using `get_cutoffs_from_model_options` and constructs the dictionary. If AtomicData_options
+        is provided, a warning is logged to indicate potential inconsistency with the model options.
+
+        The function logs the resulting AtomicData_options in a formatted and indented manner.
+
+        Parameters:
+        ----------
+            model: The model object containing model_options used to extract cutoff values.
+            AtomicData_options (dict, optional): Dictionary of atomic data options. If None, it will be generated.
+
+        Returns:
+        -------
+            dict: The updated or initialized AtomicData_options dictionary.
+        """
+        if AtomicData_options is None:
+            from dptb.utils.argcheck import get_cutoffs_from_model_options
+            # get the cutoffs from model options
+            r_max, er_max, oer_max  = get_cutoffs_from_model_options(model.model_options)
+            AtomicData_options = {'r_max': r_max, 'er_max': er_max, 'oer_max': oer_max}
+        else:
+            log.warning(msg="AtomicData_options is extracted from NEGF input file. " \
+                            "This may be not consistent with the model options. " \
+                            "Please be careful and check the cutoffs.")
+        formatted = json.dumps(AtomicData_options, indent=4)
+        indented = '\n'.join(' ' * 15 + line for line in formatted.splitlines())
+        log.info("The AtomicData_options is:\n%s", indented)
+
+        return AtomicData_options
+
 
     def fermi_dirac(self, x) -> torch.Tensor:
         return 1 / (1 + torch.exp(x / self.kBT))

dpnegf/tests/test_negf_negf_hamiltonian_init.py

@@ -75,8 +75,7 @@ def test_negf_Hamiltonian(root_directory):
             energy=e, 
             kpoint=kpoints[0], 
             eta_lead=negf_json['task_options']["eta_lead"],
-            method=negf_json['task_options']["sgf_solver"],
-            save=False
+            method=negf_json['task_options']["sgf_solver"]
             )
     print("lead_L self energy:",deviceprop.lead_L.se)
     print("lead_R self energy:",deviceprop.lead_R.se)

dpnegf/utils/argcheck.py

@@ -1047,8 +1047,8 @@ def negf():
         Argument("stru_options", dict, optional=False, sub_fields=stru_options(), doc=doc_stru_options),
         Argument("poisson_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[poisson_options()], doc=doc_poisson_options),
         Argument("sgf_solver", str, optional=True, default="Sancho-Rubio", doc=doc_sgf_solver),
-        Argument("self_energy_save", bool, optional=True, default=False, doc="whether to save the self energy"),
-        Argument("self_energy_save_path", str, optional=True, default=None, doc="the path to save the self energy"),
+        Argument("use_saved_se", bool, optional=True, default=False, doc="whether to use saved self energy"),
+        Argument("self_energy_save_path", str, optional=True, default=None, doc="the directory to save the self energy or load the self energy"),
         Argument("se_info_display", bool, optional=True, default=False, doc="whether to display the self energy information"),
         Argument("espacing", [int, float], optional=False, doc=doc_espacing),
         Argument("emin", [int, float], optional=False, doc=doc_emin),

examples/graphene/band.json

@@ -0,0 +1,18 @@
+{   
+    "structure":"train/data/POSCAR",
+    "task_options": {
+        "task": "band",
+        "kline_type":"abacus",
+        "kpath":[
+            [0, 0, 0, 50],
+            [0.5, 0, 0, 50],
+            [0.3333333, 0.3333333, 0, 50],
+            [0, 0, 0, 1]
+        ],
+        "nkpoints":151,
+        "klabels":["G", "M", "K", "G"],
+        "emin":-20,
+        "emax":20,
+        "nel_atom":{"C":4}
+    }
+}

examples/graphene/extra_baseline/c_spd.json

@@ -0,0 +1,11 @@
+{
+    "common_options": {
+        "basis": {
+                "C": [                
+                    "2s",
+                    "2p",
+                    "d*"
+                    ]
+        }
+    }
+}