Enable increasing tolerance for number of electrons in wan2skeaf.jl

* When Fermi energy doesn't converge within provided tolerance in number of electrons,
  increase the tolerance by a factor of 2 up to 0.001.
* Print a warning if the initial tolerance is higher than 0.001
* Adjust the parser accordingly, add new output parameters:
  - `tol_n_electrons_initial`
  - `tol_n_electrons_final`

Fix reference test file for SKEAF input parameters
 * Put `convert_fermi_energy_eV_to_Ry` parameter in the right place

aiida_skeaf/parsers/wan2skeaf.py

@@ -154,12 +154,15 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
     regexs = {
         "input_file_not_found": re.compile(r"ERROR: input file\s*(.+) does not exist."),
         "failed_to_find_Fermi_energy_within_tolerance": re.compile(
-            r"Error: Failed to find Fermi energy within tolerance, Δn_elec = ([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
+            r"Error: tolerance for number of electrons exceeded the tol_n_electrons_upperbound. Exiting..."
         ),
         "timestamp_started": re.compile(r"Started on\s*(.+)"),
         "num_electrons": re.compile(
             r"Number of electrons:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
         ),
+        "tol_n_electrons_initial": re.compile(
+            r"Initial tolerance for number of electrons \(default 1e-6\):\s*([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
+        ),
         "fermi_energy_in_bxsf": re.compile(
             r"Fermi Energy from file:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
         ),
@@ -185,8 +188,8 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
         "occupation_prefactor": re.compile(
             r"Occupation prefactor:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
         ),
-        "tol_n_electrons": re.compile(
-            r"Tolerance for number of electrons:\s*([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
+        "tol_n_electrons_final": re.compile(
+            r"Final tolerance for number of electrons:\s*([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
         ),
         "band_indexes_in_bxsf": re.compile(r"Bands in bxsf:\s*(.+)"),
         "timestamp_end": re.compile(r"Job done at\s*(.+)"),
@@ -232,7 +235,8 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
     ]
     float_keys = [
         "smearing_width",
-        "tol_n_electrons",
+        "tol_n_electrons_initial",
+        "tol_n_electrons_final",
         "fermi_energy_in_bxsf",
         "fermi_energy_computed",
         "fermi_energy_computed_Ry",

tests/test_data/test_input_parameters_optional.yml

@@ -1,5 +1,5 @@
-convert_fermi_energy_eV_to_Ry: false
 angle_iso_convention: false
+convert_fermi_energy_eV_to_Ry: false
 ending_phi: 90.0
 ending_theta: 90.0
 fermi_energy: 14.2832

utils/wan2skeaf.jl/wan2skeaf.jl

@@ -81,7 +81,7 @@ The script
     println("Smearing type: ", smearing_type)
     println("Smearing width: ", width_smearing)
     println("Occupation prefactor: ", prefactor)
-    println("Tolerance for number of electrons: ", tol_n_electrons)
+    println("Initial tolerance for number of electrons (default 1e-6): ", tol_n_electrons)
 
     # some times, w/o smearing, the number of electrons cannot be integrated to
     # the exact number of electrons, since we only have discrete eigenvalues.
@@ -106,14 +106,36 @@ The script
     BOHR_TO_ANG   = 0.529177210903
 
     εF_bxsf = 0.0
-    try
-        εF_bxsf = Wannier.compute_fermi_energy(eigenvalues, num_electrons, kBT, smearing; tol_n_electrons, prefactor=prefactor)
-    catch e
-        println("Error: ", e.msg)
+    tol_n_electrons_upperbound = 1e-3
+    if tol_n_electrons > tol_n_electrons_upperbound
+        println("Warning: tolerance for number of electrons is large, > 1e-3, which may lead to inaccurate Fermi energy.")
+        tol_n_electrons_upperbound = tol_n_electrons
+    end
+    tol_n_e_curr = tol_n_electrons
+
+    # try to compute Fermi energy with increasing tolerance for number of electrons
+    while tol_n_e_curr <= tol_n_electrons_upperbound
+        println("Current tolerance for number of electrons: ", tol_n_e_curr)
+        try
+            εF_bxsf = Wannier.compute_fermi_energy(eigenvalues, num_electrons, kBT, smearing; tol_n_electrons=tol_n_e_curr, prefactor=prefactor)
+            break
+        catch e
+            println("Error: ", e.msg)
+            if startswith(e.msg, "Failed to find Fermi energy within tolerance")
+                println("   Increasing tolerance for number of electrons by a factor of 2...")
+                tol_n_e_curr = tol_n_e_curr*2
+                continue
+            end
+            exit(3)
+        end
+    end
+    if tol_n_e_curr > tol_n_electrons_upperbound
+        println("Error: tolerance for number of electrons exceeded the tol_n_electrons_upperbound. Exiting...")
         exit(3)
     end
     @printf("Computed Fermi energy: %.8f\n", εF_bxsf)
     @printf("Fermi energy unit: eV\n")
+    @printf("Final tolerance for number of electrons: %.8f\n", tol_n_e_curr)
     E_bxsf = reduce(vcat, eigenvalues)
     ε_bxsf_below = maximum(E_bxsf[E_bxsf .< εF_bxsf])
     ε_bxsf_above = minimum(E_bxsf[E_bxsf .> εF_bxsf])
@@ -136,14 +158,11 @@ The script
     for ib in band_range
         # here I am still using the Fermi energy from input bxsf, i.e., QE scf Fermi
         outfile = out_filename * "_band_$(ib).bxsf"
-
-
         band_min = minimum(bxsf.E[ib:ib, :, :, :])
         band_max = maximum(bxsf.E[ib:ib, :, :, :])
-	    println("Min and max of band $ib : $band_min $band_max")
-
-	    #if (bxsf.fermi_energy >= band_min && bxsf.fermi_energy <= band_max)
-	    E_band_Ry = bxsf.E[ib:ib, :, :, :].*(ELECTRONVOLT_SI/RYDBERG_SI)
+        println("Min and max of band $ib : $band_min $band_max")
+        #if (bxsf.fermi_energy >= band_min && bxsf.fermi_energy <= band_max)
+        E_band_Ry = bxsf.E[ib:ib, :, :, :].*(ELECTRONVOLT_SI/RYDBERG_SI)
         E_fermi_Ry = bxsf.fermi_energy*(ELECTRONVOLT_SI/RYDBERG_SI)
         span_vectors_bohr = bxsf.span_vectors.*BOHR_TO_ANG/2/pi
         # what about the origin? It has to be zero (Gamma point) for bxsf so I don't change it here