New wan2skeaf in Julia with energy unit conversion from eV to Ry

aiida_skeaf/calculations/skeaf.py

@@ -145,6 +145,17 @@ def generate_input_filecontent(self) -> str:
             params["ending_theta"] = phi
             params["ending_phi"] = theta
 
+        # unit conversion, constants from QE/Modules/Constants.f90
+        ELECTRONVOLT_SI = 1.602176634e-19
+        HARTREE_SI = 4.3597447222071e-18
+        RYDBERG_SI = HARTREE_SI / 2.0
+
+        convert_fermi_energy = params.pop("convert_fermi_energy_eV_to_Ry")
+        if convert_fermi_energy:
+            params["fermi_energy"] = params["fermi_energy"] * (
+                ELECTRONVOLT_SI / RYDBERG_SI
+            )
+
         # generate the raw input for skeaf
         params = SkeafParameters(params).generate()
 

aiida_skeaf/calculations/wan2skeaf.py

@@ -5,7 +5,7 @@
 """
 import pathlib
 
-from voluptuous import Any, Optional, Required, Schema
+from voluptuous import Any, Optional, Required, Schema  # pylint: disable=unused-import
 
 from aiida import orm
 from aiida.common import datastructures
@@ -19,7 +19,7 @@ class Wan2skeafCalculation(CalcJob):
 
     _DEFAULT_INPUT_BXSF = "input.bxsf"
     _DEFAULT_OUTPUT_FILE = "wan2skeaf.out"
-    _DEFAULT_OUTPUT_BXSF = "output.bxsf"
+    _DEFAULT_OUTPUT_BXSF = "output"
 
     @classmethod
     def define(cls, spec):
@@ -105,6 +105,14 @@ def define(cls, spec):
             message="Calculation did not finish correctly.",
         )
 
+        spec.exit_code(
+            304,
+            "ERROR_NUM_ELEC_NOT_CONVERGED",
+            message="The bisection algorithm to compute Fermi level "
+            + "could not converge within the tolerance in number of electrons.\n"
+            + "Try increasing the tolerance by setting `tol_n_electrons` in the input parameters.",
+        )
+
     def prepare_for_submission(self, folder):
         """
         Create input files.
@@ -118,18 +126,22 @@ def prepare_for_submission(self, folder):
         # validate input parameters
         parameters = InputParameters(self.inputs.parameters.get_dict()).get_dict()
         cmdline_params = [
-            "--num_electrons",
+            "-n",
             parameters["num_electrons"],
-            "--band_index",
+            "-b",
             parameters["band_index"],
-            "--out_filename",
+            "-o",
             self._DEFAULT_OUTPUT_BXSF,
         ]
         if "num_spin" in parameters:
             cmdline_params += ["--num_spin", parameters["num_spin"]]
         if "smearing_type" in parameters and "smearing_value" in parameters:
-            cmdline_params += ["--smearing_type", parameters["smearing_type"]]
-            cmdline_params += ["--smearing_value", parameters["smearing_value"]]
+            cmdline_params += ["-s", parameters["smearing_type"]]
+            cmdline_params += ["-w", parameters["smearing_value"]]
+        if "occupation_prefactor" in parameters:
+            cmdline_params += ["-p", parameters["occupation_prefactor"]]
+        if "tol_n_electrons" in parameters:
+            cmdline_params += ["-t", parameters["tol_n_electrons"]]
 
         cmdline_params.append(self.inputs.bxsf_filename.value)
         #
@@ -173,9 +185,11 @@ def prepare_for_submission(self, folder):
 input_parameters = {
     Required("num_electrons"): int,
     Optional("num_spin"): int,
-    Optional("band_index", default="all"): Any(int, str),
+    Optional("band_index", default=-1): int,
     Optional("smearing_type"): str,
     Optional("smearing_value"): float,
+    Optional("occupation_prefactor"): int,
+    Optional("tol_n_electrons"): float,
 }
 
 

aiida_skeaf/data/__init__.py

@@ -29,6 +29,7 @@
     # if set as True, I will exchange theta and phi when writting the raw input file
     # for skeaf.
     Optional("angle_iso_convention", default=False): bool,
+    Optional("convert_fermi_energy_eV_to_Ry", default=False): bool,
 }
 
 # Allowed input paramters of skeaf config.in

aiida_skeaf/parsers/wan2skeaf.py

@@ -27,6 +27,15 @@ class JobNotFinishedError(
     """Raised when wan2skeaf job is not finished and end timestamp is not in the output."""
 
 
+class NumElecNotWithinToleranceError(
+    Exception
+):  # Should be inherited from aiida.common.exceptions?
+    """
+    Raised when the bisection algorithm to compute Fermi level can't converge
+    within the tolerance in number of electrons.
+    """
+
+
 class Wan2skeafParser(Parser):
     """
     Parser class for parsing output of ``wan2skeaf.py``.
@@ -73,6 +82,9 @@ def parse(self, **kwargs):
             except BXSFFileNotFoundError as exc:
                 self.logger.error(f"File not found: {exc}")
                 return self.exit_codes.ERROR_MISSING_INPUT_FILE
+            except NumElecNotWithinToleranceError as exc:
+                self.logger.error(f"Calculation failed: {exc}")
+                return self.exit_codes.ERROR_NUM_ELEC_NOT_CONVERGED
             except JobNotFinishedError as exc:
                 self.logger.error(f"Calculation not finished: {exc}")
                 return self.exit_codes.ERROR_JOB_NOT_FINISHED
@@ -96,18 +108,22 @@ def attach_bxsf_files(  # pylint: disable=inconsistent-return-statements
         """Attach RemoteData for extracted bxsf."""
 
         input_params = self.node.inputs["parameters"].get_dict()
-        input_band_index = input_params.get("band_index", "all")
+        input_band_index = input_params.get("band_index", -1)
 
-        if input_band_index == "all":
+        if input_band_index == -1:
             indexes = band_indexes_in_bxsf
         else:
             indexes = [input_band_index]
 
         remote_folder = self.node.outputs.remote_folder
         remote_folder_path = pathlib.Path(remote_folder.get_remote_path())
         remote_files = remote_folder.listdir()
-        bxsf_filename = Wan2skeafCalculation._DEFAULT_OUTPUT_BXSF.replace(  # pylint: disable=protected-access
-            ".bxsf", "_band_{:d}.bxsf"
+        # bxsf_filename = Wan2skeafCalculation._DEFAULT_OUTPUT_BXSF.replace(  # pylint: disable=protected-access
+        #     ".bxsf", "_band_{:d}.bxsf"
+        # )
+        bxsf_filename = (
+            Wan2skeafCalculation._DEFAULT_OUTPUT_BXSF
+            + "_band_{:d}.bxsf"  # pylint: disable=protected-access
         )
 
         for idx in indexes:
@@ -136,16 +152,45 @@ 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."),
+        "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: 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]+)"
         ),
         "fermi_energy_computed": re.compile(
             r"Computed Fermi energy:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
         ),
+        "fermi_energy_computed_Ry": re.compile(
+            r"Computed Fermi energy in Ry:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
+        ),
+        "fermi_energy_unit": re.compile(r"Fermi energy unit:\s*(.+)"),
+        "closest_eigenvalue_below_fermi": re.compile(
+            r"Closest eigenvalue below Fermi energy:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
+        ),
+        "closest_eigenvalue_above_fermi": re.compile(
+            r"Closest eigenvalue above Fermi energy:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
+        ),
         "num_bands": re.compile(r"Number of bands:\s*([0-9]+)"),
         "kpoint_mesh": re.compile(r"Grid shape:\s*(.+)"),
+        "smearing_type": re.compile(r"Smearing type:\s*(.+)"),
+        "smearing_width": re.compile(
+            r"Smearing width:\s*([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
+        ),
+        "occupation_prefactor": re.compile(
+            r"Occupation prefactor:\s*([+-]?(?:[0-9]*[.])?[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*(.+)"),
     }
@@ -176,15 +221,36 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
         raise BXSFFileNotFoundError(
             errno.ENOENT, os.strerror(errno.ENOENT), parameters["input_file_not_found"]
         )
+    if "failed_to_find_Fermi_energy_within_tolerance" in parameters:
+        raise NumElecNotWithinToleranceError(
+            "Failed to find Fermi energy within tolerance, Δn_elec = "
+            + f"{parameters['failed_to_find_Fermi_energy_within_tolerance']}"
+        )
     if "timestamp_end" not in parameters:
         raise JobNotFinishedError("Job not finished!")
 
     parameters["kpoint_mesh"] = [int(_) for _ in parameters["kpoint_mesh"].split("x")]
     parameters["band_indexes_in_bxsf"] = [
         int(_) for _ in parameters["band_indexes_in_bxsf"].split()
     ]
-    parameters["fermi_energy_in_bxsf"] = float(parameters["fermi_energy_in_bxsf"])
-    parameters["fermi_energy_computed"] = float(parameters["fermi_energy_computed"])
+    float_keys = [
+        "smearing_width",
+        "tol_n_electrons_initial",
+        "tol_n_electrons_final",
+        "fermi_energy_in_bxsf",
+        "fermi_energy_computed",
+        "fermi_energy_computed_Ry",
+        "closest_eigenvalue_below_fermi",
+        "closest_eigenvalue_above_fermi",
+    ]
+    for key in float_keys:
+        parameters[key] = float(parameters[key])
+    parameters["fermi_energy_unit"] = parameters["fermi_energy_unit"]
+    parameters["smearing_type"] = parameters["smearing_type"]
+    parameters["num_bands"] = int(parameters["num_bands"])
+    parameters["num_electrons"] = int(parameters["num_electrons"])
+    parameters["occupation_prefactor"] = int(parameters["occupation_prefactor"])
+
     # make sure the order is the same as parameters["band_indexes_in_bxsf"]
     parameters["band_min"] = [
         band_minmax[_][0] for _ in parameters["band_indexes_in_bxsf"]

aiida_skeaf/workflows/protocols/skeaf.yaml

@@ -10,7 +10,7 @@ default_inputs:
                 withmpi: False
         parameters:
             # num_electrons: 1
-            band_index: 'all'
+            band_index: -1
     skeaf:
         metadata:
             options:
@@ -35,6 +35,7 @@ default_inputs:
             ending_phi: 0.000000
             num_rotation: 45
             angle_iso_convention: True
+            convert_fermi_energy_eV_to_Ry: True
 default_protocol: moderate
 protocols:
     moderate:

examples/example_02.py

@@ -9,20 +9,23 @@
 def submit():
     """Submit a Wan2skeafCalculation."""
 
-    bxsf = orm.load_node(137298)
+    # bxsf = orm.load_node(137298)
+    bxsf = orm.load_node("e1f4857d-e6d0-4a7f-a25c-783d2b9ca24a")
 
     # skeaf_code = orm.load_code("skeaf@localhost")
     wan2skeaf_code = orm.load_code("skeaf-wan2skeaf@localhost")
 
     parameters = {
-        "num_electrons": 1,
-        "band_index": "all",
+        "num_electrons": 11,
+        "band_index": -1,
     }
 
     inputs = {
         "code": wan2skeaf_code,
         "parameters": parameters,
         "bxsf": bxsf,
+        "bxsf_filename": "bxsf.7z",
+        "settings": dict({"autolink_bxsf_filename": "bxsf.7z"}),
     }
 
     calc = engine.submit(Wan2skeafCalculation, **inputs)

examples/example_03.py

@@ -4,35 +4,43 @@
 
 from aiida_wannier90_workflows.utils.workflows.builder.serializer import print_builder
 
-from aiida_skeaf.utils.plot import plot_frequency_workchain
+from aiida_skeaf.utils.plot import (  # pylint: disable=unused-import
+    plot_frequency_workchain,
+)
 from aiida_skeaf.workflows import SkeafWorkChain
 
 
 def submit():
     """Submit a SkeafWorkChain."""
 
-    bxsf = orm.load_node(137298)
+    # bxsf = orm.load_node(137298)
+    bxsf = orm.load_node("e1f4857d-e6d0-4a7f-a25c-783d2b9ca24a")
 
     codes = {
-        "wan2skeaf": "skeaf-wan2skeaf@localhost",
-        "skeaf": "skeaf@localhost",
+        "wan2skeaf": "skeaf-wan2skeaf-jl@localhost",
+        "skeaf": "skeaf-skeaf@localhost",
     }
 
     builder = SkeafWorkChain.get_builder_from_protocol(
         codes,
         bxsf=bxsf,
-        num_electrons=2,
+        num_electrons=11,
         protocol="fast",
     )
 
+    builder.wan2skeaf.settings.autolink_bxsf_filename = (  # pylint: disable=no-member
+        "bxsf.7z"
+    )
+    builder.wan2skeaf.bxsf_filename = "bxsf.7z"  # pylint: disable=no-member
+
     print_builder(builder)
 
     wkchain = engine.submit(builder)
 
     print(f"Submitted {wkchain}")
 
     # Once workchain has finished, plot frequency vs angle
-    plot_frequency_workchain(wkchain, x="theta")
+    # plot_frequency_workchain(wkchain, x="theta")
 
 
 if __name__ == "__main__":

examples/example_05.py

@@ -1,32 +1,39 @@
 #!/usr/bin/env runaiida
 """Submit a wan2skeaf calculation with cold smearing."""
 # pylint: disable=unused-import
+import pathlib
 
 from aiida import engine, orm
 
+from aiida_skeaf.calculations.functions import create_bxsf_from_file
 from aiida_skeaf.calculations.wan2skeaf import Wan2skeafCalculation
 
+INPUT_DIR = pathlib.Path(__file__).absolute().parent / "input_files"
+
 
 def submit():
     """Submit a Wan2skeafCalculation."""
 
-    bxsf = orm.load_node(152164)
+    remote_path = orm.Str(str(INPUT_DIR / "bxsf"))
+    computer = orm.Str("localhost")
+    bxsf = create_bxsf_from_file(remote_path, computer)
+    # bxsf = orm.load_node("e1f4857d-e6d0-4a7f-a25c-783d2b9ca24a")
 
-    wan2skeaf_code = orm.load_code("skeaf-wan2skeaf@localhost-slurm")
+    wan2skeaf_code = orm.load_code("skeaf-wan2skeaf-jl@localhost")
 
     parameters = {
-        "num_electrons": 31,
-        "band_index": "all",
+        "num_electrons": 11,
+        "band_index": -1,
         "smearing_type": "cold",
-        "smearing_value": 0.13605693122994,
+        "smearing_value": 0.1,
     }
 
     inputs = {
         "code": wan2skeaf_code,
         "parameters": parameters,
         "bxsf": bxsf,
-        "bxsf_filename": "nscf.bxsf",
-        "settings": dict({"autolink_bxsf_filename": "nscf.bxsf"}),
+        "bxsf_filename": "copper.7z",
+        "settings": dict({"autolink_bxsf_filename": "copper.7z"}),
     }
 
     calc = engine.submit(Wan2skeafCalculation, **inputs)