Refactor: move example files into examples/ dir and add validation

SIAB/io/read_input.py

@@ -13,7 +13,7 @@ def read_siab_plaintext(fname: str = ""):
     float_pattern = r"^\d+\.\d*$"
     int_pattern = r"^\d+$"
     scalar_keywords = ["Ecut", "sigma", "element"]
-    result = {}
+    result = {"fit_basis": "pw"}
     if fname == "":
         raise ValueError("No filename provided")
     with open(fname, "r") as f:
@@ -471,6 +471,11 @@ def skip_ppread(user_settings: dict):
     """
 
     skip = True
+    # case 0
+    # if element is not set, it is not possible to skip the pseudopotential read-in
+    if "element" not in user_settings:
+        print("AUTOSET: `element` is not specified => AUTOSET", flush=True)
+        return False
     # case 1
     # if nbands is specified as auto, occ, all, it must requires the number of valence
     # electrons, therefore it is not possible to skip the pseudopotential read-in
@@ -536,10 +541,51 @@ def skip_ppread(user_settings: dict):
         # in a list in input, therefore the index is known).
     return skip
 
+def _validate_param(user_settings: dict):
+    """validate the input parameters
+    
+    Parameters
+    ----------
+    user_settings: dict
+        the user settings
+    
+    Returns
+    -------
+    None
+    """
+    # check if the shape assigned to orbitals is valid
+    shape2index = {rs["shape"]: i for i, rs in enumerate(user_settings["reference_systems"])}
+    for iorb, orb in enumerate(user_settings["orbitals"]):
+        shape = orb["shape"]
+        shape = [shape] if not isinstance(shape, list) else shape
+        for s in shape:
+            assert isinstance(s, (str, int)), f"shape {s} is not a valid shape"
+            if isinstance(s, str):
+                assert s in shape2index, f"shape {s} is not found in reference systems"
+
+    # check if the nbands set for reference system is smaller than
+    # bands needed for fitting orbitals
+    shape2index = {rs["shape"]: i for i, rs in enumerate(user_settings["reference_systems"])}
+    for iorb, orb in enumerate(user_settings["orbitals"]):
+        shape = orb["shape"]
+        shape = [shape] if not isinstance(shape, list) else shape
+        for s in shape:
+            if isinstance(s, str):
+                assert orb["nbands_ref"] <= user_settings["reference_systems"][shape2index[s]]["nbands"], \
+                    f"ERROR: `nbands_ref` for orbital {iorb} is larger than the number of bands set for\
+ reference system `{s}`"
+            elif isinstance(s, int):
+                assert orb["nbands_ref"] <= user_settings["reference_systems"][s]["nbands"], \
+                    f"ERROR: `nbands_ref` for orbital {iorb} is larger than the number of bands set for\
+ reference system `{s}`"
+
+
 def parse(user_settings: dict):
     """unpack the SIAB input to structure (shape as key and bond lengths are list as value),
     input setting of abacus, orbital generation settings, environmental settings and general description
     """
+    _validate_param(user_settings)
+
     # move the information fetch from pseudopotential from front.py here...
     # get value from the dict returned by function from_pseudopotential
 

SIAB/spillage/api.py

@@ -686,6 +686,9 @@ def _nzeta_infer(folder, nband):
     for isk in range(nspin*len(wk)): # loop over (ispin, ik)
         w = wk[isk % len(wk)] # spin-up and spin-down share the wk
         wfc, _, _, _ = read_wfc_lcao_txt(os.path.join(outdir, f"{fwfc}{isk+1}.txt"))
+        assert wfc.shape[1] >= nband, \
+            f"ERROR: number of bands for orbgen is larger than calculated: {nband} > {wfc.shape[1]}"
+
         # the complete return list is (wfc.T, e, occ, k)
         ovlp = read_triu(os.path.join(outdir, f"data-{isk}-S"))
 

examples/jy.json

@@ -0,0 +1,56 @@
+{
+    "environment": "",
+    "mpi_command": "mpirun -np 8",
+    "abacus_command": "abacus",
+
+    "pseudo_dir": "/root/abacus-develop/pseudopotentials/sg15_oncv_upf_2020-02-06/",
+    "pseudo_name": "Si_ONCV_PBE-1.0.upf",
+    "ecutwfc": 100,
+    "bessel_nao_smooth": 0,
+    "bessel_nao_rcut": [6, 7, 8, 9, 10],
+    "smearing_sigma": 0.01,
+
+    "fit_basis": "jy",
+    "optimizer": "bfgs",
+    "max_steps": 3000,
+    "spill_guess": "atomic",
+    "nthreads_rcut": 4,
+    "jY_type": "reduced",
+
+    "reference_systems": [
+        {
+            "shape": "dimer",
+            "nbands": 8,
+            "nspin": 1,
+            "lmaxmax": 2,
+            "bond_lengths": [1.62, 1.82, 2.22, 2.72, 3.22]
+        },
+        {
+            "shape": "trimer",
+            "nbands": 10,
+            "nspin": 1,
+            "lmaxmax": 2,
+            "bond_lengths": [1.9, 2.1, 2.6]
+        }
+    ],
+    "orbitals": [
+        {
+            "zeta_notation": "auto",
+            "shape": "dimer",
+            "nbands_ref": 4,
+            "orb_ref": "none"
+        },
+        {
+            "zeta_notation": "auto",
+            "shape": "dimer",
+            "nbands_ref": 10,
+            "orb_ref": 0
+        },
+        {
+            "zeta_notation": "auto",
+            "shape": "trimer",
+            "nbands_ref": 20,
+            "orb_ref": 1
+        }
+    ]
+}

SIAB_INPUT.json → examples/pw.json

@@ -12,7 +12,7 @@
     "pseudo_rcut": 10,
     "pseudo_mesh": 1,
 
-    "basis_type": "jy",
+    "fit_basis": "pw",
 
     "optimizer": "bfgs",
     "max_steps": 1000,