Skip to content

Commit

Permalink
Catch the error when Fermi energy calculation couldn't converge withi…
Browse files Browse the repository at this point in the history 
…n the tolerance in number of electrons

- Catch the error in `wan2skeaf.jl` and print it to the output file
- Parse the error message and add new exit code (304) to Wan2skeafCalculation
- Pass `tol_n_elec` parameter to `wan2skeaf.jl` suth that the tolerance in number of electrons can be controlled
- Make `wan2skeaf.jl` output more verbose, print also the parameters that are actually used in the calculation
- Parse additional wan2skeaf output to the `output_parameters`
  • Loading branch information
@npaulish @qiaojunfeng
npaulish authored and qiaojunfeng committed Jun 4, 2024
1 parent af18047 commit 31962c0
Show file tree
Hide file tree
Showing 3 changed files with 90 additions and 24 deletions.
11 changes: 11 additions & 0 deletions aiida_skeaf/calculations/wan2skeaf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -132,6 +140,8 @@ def prepare_for_submission(self, folder):
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)
#
Expand Down Expand Up @@ -179,6 +189,7 @@ def prepare_for_submission(self, folder):
Optional("smearing_type"): str,
Optional("smearing_value"): float,
Optional("occupation_prefactor"): int,
Optional("tol_n_electrons"): float,
}


Expand Down
62 changes: 50 additions & 12 deletions aiida_skeaf/parsers/wan2skeaf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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``.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -140,8 +152,14 @@ 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: Failed to find Fermi energy within tolerance, Δn_elec = ([+-]?(?:[0-9]*[.])?[0-9]+e?[+-]?[0-9]*)"
),
"timestamp_started": re.compile(r"Started on\s*(.+)"),
"num_electrons": re.compile(
r"Number of electrons:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
),
"fermi_energy_in_bxsf": re.compile(
r"Fermi Energy from file:\s*([+-]?(?:[0-9]*[.])?[0-9]+)"
),
Expand All @@ -160,6 +178,16 @@ def parse_wan2skeaf_out(filecontent: ty.List[str]) -> orm.Dict:
),
"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": re.compile(
r"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*(.+)"),
}
Expand Down Expand Up @@ -190,25 +218,35 @@ 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"])
parameters["fermi_energy_computed_Ry"] = float(
parameters["fermi_energy_computed_Ry"]
)
float_keys = [
"smearing_width",
"tol_n_electrons",
"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["closest_eigenvalue_below_fermi"] = float(
parameters["closest_eigenvalue_below_fermi"]
)
parameters["closest_eigenvalue_above_fermi"] = float(
parameters["closest_eigenvalue_above_fermi"]
)
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"]
Expand Down
41 changes: 29 additions & 12 deletions utils/wan2skeaf.jl/wan2skeaf.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,30 +38,35 @@ The script
- `-s, --smearing_type`: smearing type, default is `NoneSmearing()` (no smearing)
- `-w, --width_smearing`: smearing width, default is 0.0
- `-p, --prefactor`: occupation prefactor, 2 for non SOC, 1 for SOC, default is 2
- `-t, --tol_n_electrons`: tolerance for number of electrons, default is 1e-6
"""
@main function main(bxsf::String; num_electrons::Int, band_index::Int=-1, out_filename::String="skeaf", smearing_type::String="none", width_smearing::Float64=0.0, prefactor::Int=2)
@main function main(bxsf::String; num_electrons::Int, band_index::Int=-1, out_filename::String="skeaf", smearing_type::String="none", width_smearing::Float64=0.0, prefactor::Int=2, tol_n_electrons::Float64=1e-6)
println("Started on ", Dates.now())
if !isfile(bxsf)
println("ERROR: Input file $bxsf does not exist.")
println("ERROR: input file $bxsf does not exist.")
exit(2)
end
if endswith(bxsf, ".7z")
# -y: assume yes, will overwrite existing unzipped files
ret = run(`$(p7zip()) x -y $bxsf`)
success(ret) || error("Error while unzipping file: $bxsf")
success(ret) || error("error while unzipping file: $bxsf")
bxsf_files = filter(x -> endswith(x, ".bxsf"), readdir("."))
if length(bxsf_files) != 1
error("Expecting only one .bxsf file in the 7z file, but got $(length(bxsf_files))")
error("expecting only one .bxsf file in the 7z file, but got $(length(bxsf_files))")
end
bxsf_filename = bxsf_files[1]
dst_filename = "input.bxsf"
mv(bxsf_filename, dst_filename) # not sure this is needed! skeaf reads the output file?
bxsf = dst_filename
end
println("Number of electrons: ", num_electrons)
bxsf = WannierIO.read_bxsf(bxsf)
@printf("Fermi Energy from file: %.8f\n", bxsf.fermi_energy)

nbands, nkx, nky, nkz = size(bxsf.E)
println("Number of bands: ", nbands)
println("Grid shape: $nkx x $nky x $nkz")

kBT = width_smearing
if smearing_type == "none"
smearing = Wannier.NoneSmearing()
Expand All @@ -70,12 +75,13 @@ The script
elseif smearing_type == "marzari-vanderbilt" || smearing_type == "cold"
smearing = Wannier.ColdSmearing()
else
error("Unknown smearing type: $smearing_type")
error("unknown smearing type: $smearing_type")
end

nbands, nkx, nky, nkz = size(bxsf.E)
println("Number of bands: ", nbands)
println("Grid shape: $nkx x $nky x $nkz")
println("Smearing type: ", smearing_type)
println("Smearing width: ", width_smearing)
println("Occupation prefactor: ", prefactor)
println("Tolerance for number of electrons: ", 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.
Expand All @@ -84,7 +90,7 @@ The script
# always possible to integrate to the exact integer for number of electrons.
# tol_n_electrons = 1e-3
# this is the default
tol_n_electrons = 1e-6
# tol_n_electrons = 1e-6

# note that according to bxsf specification, the eigenvalues are defined
# on "general grid" instead of "periodic grid", i.e., the right BZ edges are
Expand All @@ -99,16 +105,27 @@ The script
RYDBERG_SI = HARTREE_SI/2.0
BOHR_TO_ANG = 0.529177210903

εF_bxsf = Wannier.compute_fermi_energy(eigenvalues, num_electrons, kBT, smearing; tol_n_electrons, prefactor=prefactor)
ε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)
exit(3)
end
@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)
@printf("Closest eigenvalue above Fermi energy: %.8f\n", ε_bxsf_above)

@printf("Computed Fermi energy in Ry: %.8f\n", εF_bxsf*(ELECTRONVOLT_SI/RYDBERG_SI))
println("Constants used for the conversion (from QE/Modules/Constants.f90): ")
println(" ELECTRONVOLT_SI: ", ELECTRONVOLT_SI)
println(" RYDBERG_SI: ", RYDBERG_SI)
println(" BOHR_TO_ANG: ", BOHR_TO_ANG)

# write each band into one bxsf file
if band_index < 0
band_range = 1:nbands
Expand Down

0 comments on commit 31962c0

Please sign in to comment.