Conversion from eV to Ry only for SKEAF input files.

- Modify wan2skeaf.jl such that only the energies written to bxsfs for SKEAF inputs are converted from eV to Ry.
- Fermi energy and min and max band energies are kept in eV in the output for convenience (the units of the original bxsf file).
- Fermi enery in Ry is also printed to the output.
- Add `convert_fermi_energy_eV_to_Ry` flag to the input parameters for the SkeafCalculation and set it to `True` for default protocols.
- This will allow directly specifying Fermi energy in eV when setting up a SkeafWorkChain, then the conversion will be done when preparing the SKEAF input file.
- If Fermi energy for SKEAF calculations is specified in Ry, `convert_fermi_energy_eV_to_Ry` should be set to `False` to skip conversion.

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/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

@@ -148,8 +148,8 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
         "fermi_energy_computed": re.compile(
             r"Computed Fermi energy:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
         ),
-        "fermi_energy_computed_eV": re.compile(
-            r"Computed Fermi energy in eV:\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(
@@ -199,8 +199,8 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
     ]
     parameters["fermi_energy_in_bxsf"] = float(parameters["fermi_energy_in_bxsf"])
     parameters["fermi_energy_computed"] = float(parameters["fermi_energy_computed"])
-    parameters["fermi_energy_computed_eV"] = float(
-        parameters["fermi_energy_computed_eV"]
+    parameters["fermi_energy_computed_Ry"] = float(
+        parameters["fermi_energy_computed_Ry"]
     )
     parameters["fermi_energy_unit"] = parameters["fermi_energy_unit"]
     parameters["closest_eigenvalue_below_fermi"] = float(

aiida_skeaf/workflows/protocols/skeaf.yaml

@@ -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:

tests/test_data/test_input_parameters_optional.yml

@@ -1,3 +1,4 @@
+convert_fermi_energy_eV_to_Ry: false
 angle_iso_convention: false
 ending_phi: 90.0
 ending_theta: 90.0

utils/wan2skeaf.jl/wan2skeaf.jl

@@ -100,14 +100,14 @@ The script
     BOHR_TO_ANG   = 0.529177210903
 
     εF_bxsf = Wannier.compute_fermi_energy(eigenvalues, num_electrons, kBT, smearing; tol_n_electrons, prefactor=prefactor)
-    @printf("Computed Fermi energy: %.8f\n", εF_bxsf*(ELECTRONVOLT_SI/RYDBERG_SI))
-    @printf("Computed Fermi energy in eV: %.8f\n", εF_bxsf)
-    @printf("Fermi energy unit: Ry\n")
+    @printf("Computed Fermi energy: %.8f\n", εF_bxsf)
+    @printf("Computed Fermi energy in Ry: %.8f\n", εF_bxsf*(ELECTRONVOLT_SI/RYDBERG_SI))
+    @printf("Fermi energy unit: eV\n")
     E_bxsf = reduce(vcat, eigenvalues)
     ε_bxsf_below = maximum(E_bxsf[E_bxsf .< εF_bxsf])
     ε_bxsf_above = minimum(E_bxsf[E_bxsf .> εF_bxsf])
-    @printf("Closest eigenvalue below Fermi energy: %.8f\n", ε_bxsf_below *(ELECTRONVOLT_SI/RYDBERG_SI))
-    @printf("Closest eigenvalue above Fermi energy: %.8f\n", ε_bxsf_above *(ELECTRONVOLT_SI/RYDBERG_SI))
+    @printf("Closest eigenvalue below Fermi energy: %.8f\n", ε_bxsf_below)
+    @printf("Closest eigenvalue above Fermi energy: %.8f\n", ε_bxsf_above)
 
     # write each band into one bxsf file
     if band_index < 0
@@ -121,8 +121,8 @@ The script
         outfile = out_filename * "_band_$(ib).bxsf"
 
 
-        band_min = minimum(bxsf.E[ib:ib, :, :, :])*(ELECTRONVOLT_SI/RYDBERG_SI)
-        band_max = maximum(bxsf.E[ib:ib, :, :, :])*(ELECTRONVOLT_SI/RYDBERG_SI)
+        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)