diff --git a/src/aiida_wannier90_workflows/utils/workflows/plot/bands.py b/src/aiida_wannier90_workflows/utils/workflows/plot/bands.py index 2c12a5c6..17e6faaa 100755 --- a/src/aiida_wannier90_workflows/utils/workflows/plot/bands.py +++ b/src/aiida_wannier90_workflows/utils/workflows/plot/bands.py @@ -379,6 +379,17 @@ def get_workchain_fermi_energy( else: raise ValueError("Cannot find fermi energy") + if fermi_energy is None: + # `get_fermi_energy_from_nscf` (used by the nscf branches above) can + # return None. Fail loudly here rather than letting None flow into the + # band-plotting arithmetic downstream. + raise ValueError( + f"Could not read a Fermi energy from {workchain}. Its " + "`output_parameters` must contain `fermi_energy` (or both " + "`fermi_energy_up` and `fermi_energy_down`), together with " + "`fermi_energy_units` set to `eV`." + ) + return fermi_energy diff --git a/src/aiida_wannier90_workflows/utils/workflows/pw.py b/src/aiida_wannier90_workflows/utils/workflows/pw.py index f4b3635c..b6fe4870 100644 --- a/src/aiida_wannier90_workflows/utils/workflows/pw.py +++ b/src/aiida_wannier90_workflows/utils/workflows/pw.py @@ -28,13 +28,20 @@ def get_fermi_energy(output_parameters: orm.Dict) -> ty.Optional[float]: def get_fermi_energy_from_nscf( calc_nscf: ty.Union[PwBaseWorkChain, PwCalculation] -) -> float: - """Parse nscf output to get the scf Fermi energy. +) -> ty.Optional[float]: + """Get a Fermi energy from an nscf run. + + Prefer the scf Fermi energy reported in the nscf stdout via the + ``(compare with: ... computed in scf)`` marker. When that marker is absent + (see the fallback below), return the nscf's own Fermi energy instead, taken + from the parsed ``output_parameters``. The two can differ; the scf value is + kept as the first choice to preserve existing behaviour. :param calc_nscf: a nscf PwBaseWorkChain or PwCalculation :type calc_nscf: ty.Union[PwBaseWorkChain, PwCalculation] - :return: scf Fermi energy - :rtype: float + :return: the scf Fermi energy if the stdout marker is present, otherwise the + nscf Fermi energy from ``output_parameters``, else None. Unit is eV. + :rtype: float, None """ import re @@ -69,4 +76,28 @@ def get_fermi_energy_from_nscf( fermi_energy = float(match.group(1)) break + if fermi_energy is None: + # The regex above only matches the single-value "(compare with: X eV, + # computed in scf)" marker. QE prints that marker only for a metallic + # (smearing/tetrahedra), single-Fermi-energy nscf run: insulators print + # HOMO/LUMO instead, and constrained-magnetization runs print a + # two-value variant the regex does not match (see QE + # PW/src/print_ks_energies.f90). In those cases fall back to the Fermi + # energy the parser stored from this nscf run. + output_parameters = calc_nscf.outputs.output_parameters.get_dict() + # The aiida-quantumespresso parser always stores Fermi energies in eV, + # but guard on the units regardless, to mirror `get_fermi_energy` and + # avoid silently returning a value in the wrong unit. + if output_parameters.get("fermi_energy_units") == "eV": + fermi_energy = output_parameters.get("fermi_energy") + if fermi_energy is None: + # Spin-polarised runs with a constrained total magnetization + # report one Fermi level per channel and have no single chemical + # potential. Take the higher of the two as a conservative + # reference for the (frozen) energy windows. + up = output_parameters.get("fermi_energy_up") + down = output_parameters.get("fermi_energy_down") + if up is not None and down is not None: + fermi_energy = max(up, down) + return fermi_energy diff --git a/src/aiida_wannier90_workflows/workflows/wannier90.py b/src/aiida_wannier90_workflows/workflows/wannier90.py index c35e5bdb..3d00440a 100644 --- a/src/aiida_wannier90_workflows/workflows/wannier90.py +++ b/src/aiida_wannier90_workflows/workflows/wannier90.py @@ -826,20 +826,33 @@ def prepare_wannier90_pp_inputs(self): # pylint: disable=too-many-statements # Add Fermi energy if "workchain_scf" in self.ctx: - scf_output_parameters = self.ctx.workchain_scf.outputs.output_parameters - fermi_energy = get_fermi_energy(scf_output_parameters) - elif "workchain_nscf" in self.ctx: - if ( - "fermi_energy" not in parameters - ): # we can provide it if the workchain_nscf was already performed before this run. - fermi_energy = get_fermi_energy_from_nscf(self.ctx.workchain_nscf) - else: - fermi_energy = parameters["fermi_energy"] + fermi_source = self.ctx.workchain_scf + fermi_energy = get_fermi_energy(fermi_source.outputs.output_parameters) + elif "workchain_nscf" in self.ctx and "fermi_energy" not in parameters: + fermi_source = self.ctx.workchain_nscf + fermi_energy = get_fermi_energy_from_nscf(fermi_source) + elif "fermi_energy" in parameters: + # Given by the caller, which is the only source when the scf and nscf + # ran before this workchain. + fermi_source = None + fermi_energy = parameters["fermi_energy"] else: - if "fermi_energy" in parameters: - fermi_energy = parameters["fermi_energy"] - else: - raise ValueError("Cannot retrieve Fermi energy from scf or nscf output") + raise ValueError("Cannot retrieve Fermi energy from scf or nscf output") + if fermi_energy is None: + # Fail loudly here rather than passing None through to the .win + # writer, which rejects it with an opaque "Invalid value" error. + if fermi_source is None: + raise ValueError( + "The `fermi_energy` in the wannier90 parameters is None. Set it " + "to a number in eV, or drop it to read it from the scf or nscf " + "output." + ) + raise ValueError( + f"Could not read a Fermi energy from {fermi_source}. Its " + "`output_parameters` must contain `fermi_energy` (or both " + "`fermi_energy_up` and `fermi_energy_down`), together with " + "`fermi_energy_units` set to `eV`." + ) parameters["fermi_energy"] = fermi_energy inputs.parameters = orm.Dict(parameters) diff --git a/tests/utils/workflows/test_pw.py b/tests/utils/workflows/test_pw.py new file mode 100644 index 00000000..e63b8436 --- /dev/null +++ b/tests/utils/workflows/test_pw.py @@ -0,0 +1,116 @@ +"""Unit tests for :py:mod:`aiida_wannier90_workflows.utils.workflows.pw`.""" + +import pytest + +from aiida_quantumespresso.calculations.pw import PwCalculation + +from aiida_wannier90_workflows.utils.workflows.pw import get_fermi_energy_from_nscf + +# stdout of an nscf run that still prints the scf-Fermi companion marker. +NSCF_STDOUT_WITH_MARKER = """ + End of band structure calculation + + the Fermi energy is 5.9816 ev + (compare with: 5.9034 eV, computed in scf) + + Writing output data file aiida.save +""" + +# stdout of an nscf run that only prints its own Fermi energy (no marker). +NSCF_STDOUT_WITHOUT_MARKER = """ + End of band structure calculation + + the Fermi energy is 5.9816 ev + + Writing output data file aiida.save +""" + + +class _FakeRetrieved: + """Stand in for the ``retrieved`` FolderData node.""" + + def __init__(self, content): + self._content = content + + def get_object_content(self, name): # pylint: disable=unused-argument + return self._content + + +class _FakeDict: + """Stand in for an ``orm.Dict`` output node.""" + + def __init__(self, dictionary): + self._dictionary = dictionary + + def get_dict(self): + return dict(self._dictionary) + + +class _FakeOutputs: + def __init__(self, stdout, output_parameters): + self.retrieved = _FakeRetrieved(stdout) + self.output_parameters = _FakeDict(output_parameters) + + +class _FakeNscfCalc: + """Minimal stub of a finished nscf ``PwCalculation`` node. + + Only the attributes accessed by ``get_fermi_energy_from_nscf`` are + implemented, so the test needs no AiiDA profile or database. + """ + + process_class = PwCalculation + is_finished_ok = True + + def __init__(self, stdout, output_parameters): + self.outputs = _FakeOutputs(stdout, output_parameters) + + +def test_get_fermi_energy_from_nscf_marker_present(): + """Marker present: the scf value from stdout is used, not the fallback.""" + calc = _FakeNscfCalc( + NSCF_STDOUT_WITH_MARKER, + # A different value here would be returned only if the fallback ran. + {"fermi_energy": 7.0, "fermi_energy_units": "eV"}, + ) + assert get_fermi_energy_from_nscf(calc) == pytest.approx(5.9034) + + +def test_get_fermi_energy_from_nscf_fallback_to_parsed_value(): + """Marker absent: fall back to the parsed nscf Fermi energy.""" + calc = _FakeNscfCalc( + NSCF_STDOUT_WITHOUT_MARKER, + {"fermi_energy": 5.9816, "fermi_energy_units": "eV"}, + ) + assert get_fermi_energy_from_nscf(calc) == pytest.approx(5.9816) + + +def test_get_fermi_energy_from_nscf_fallback_spin_polarised(): + """Marker absent, two Fermi levels: return the higher of the two.""" + calc = _FakeNscfCalc( + NSCF_STDOUT_WITHOUT_MARKER, + { + "fermi_energy_up": 5.5, + "fermi_energy_down": 6.1, + "fermi_energy_units": "eV", + }, + ) + assert get_fermi_energy_from_nscf(calc) == pytest.approx(6.1) + + +def test_get_fermi_energy_from_nscf_returns_none_when_unavailable(): + """Marker absent and nothing parsed: return None (guarded by callers).""" + calc = _FakeNscfCalc( + NSCF_STDOUT_WITHOUT_MARKER, + {"fermi_energy_units": "eV"}, + ) + assert get_fermi_energy_from_nscf(calc) is None + + +def test_get_fermi_energy_from_nscf_ignores_non_ev_units(): + """A parsed Fermi energy in non-eV units is not used by the fallback.""" + calc = _FakeNscfCalc( + NSCF_STDOUT_WITHOUT_MARKER, + {"fermi_energy": 0.44, "fermi_energy_units": "Ry"}, + ) + assert get_fermi_energy_from_nscf(calc) is None diff --git a/tests/workflows/test_wannier90.py b/tests/workflows/test_wannier90.py index bc1a5568..cf0a2b5c 100644 --- a/tests/workflows/test_wannier90.py +++ b/tests/workflows/test_wannier90.py @@ -3,12 +3,26 @@ import io from plumpy.process_states import ProcessState +import pytest from aiida import orm from aiida.common import LinkType from aiida_quantumespresso.calculations.helpers import pw_input_helper +# The scf Fermi energy an nscf run reports in its stdout. Any test that reads a +# different value did not consult the nscf. +NSCF_SCF_FERMI_ENERGY = 5.9034 + +NSCF_STDOUT = f""" + End of band structure calculation + + the Fermi energy is 5.9816 ev + (compare with: {NSCF_SCF_FERMI_ENERGY} eV, computed in scf) + + Writing output data file aiida.save +""" + def test_scdm( generate_workchain_wannier90, @@ -35,7 +49,11 @@ def test_scdm( scf_workchain, link_type=LinkType.RETURN, link_label="remote_folder" ) - params = orm.Dict({"fermi_energy": 6.0, "number_of_electrons": 8}) + # The pw.x parser always stores `fermi_energy_units` alongside + # `fermi_energy`; `get_fermi_energy` reads a value only when it does. + params = orm.Dict( + {"fermi_energy": 6.0, "fermi_energy_units": "eV", "number_of_electrons": 8} + ) params.store() params.base.links.add_incoming( scf_workchain, link_type=LinkType.RETURN, link_label="output_parameters" @@ -108,6 +126,9 @@ def test_scdm( # mock run wannier90 pp w90pp_workchain = workchain.run_wannier90_pp()["workchain_wannier90_pp"] + # The scf Fermi energy reaches the wannier90 parameters + assert w90pp_workchain.inputs.wannier90.parameters["fermi_energy"] == 6.0 + # The wannier90 step will use `get_last_calcjob` to retrieve input parameters of the calcjob entry_point_calc_job = "wannier90.wannier90" calcjob = generate_calc_job_node( @@ -191,3 +212,110 @@ def test_scdm( _ in workchain.outputs for _ in ("scf", "nscf", "projwfc", "wannier90_pp", "pw2wannier90", "wannier90") ) + + +def _generate_finished_nscf_calc(generate_calc_job_node, fixture_localhost): + """Return a finished nscf ``PwCalculation`` node whose stdout reports the scf Fermi energy.""" + node = generate_calc_job_node("quantumespresso.pw", fixture_localhost, store=False) + node.set_process_state(ProcessState.FINISHED) + node.set_exit_status(0) + node.store() + + retrieved = orm.FolderData() + retrieved.put_object_from_filelike(io.StringIO(NSCF_STDOUT), "aiida.out") + retrieved.base.links.add_incoming( + node, link_type=LinkType.CREATE, link_label="retrieved" + ) + retrieved.store() + + return node + + +def _generate_workchain_without_scf_context( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, + w90_parameters, +): + """Return a ``Wannier90WorkChain`` whose context holds an nscf run but no scf run. + + ``w90_parameters`` are the wannier90 parameters the caller supplies, which is + where an externally computed Fermi energy would arrive. + """ + inputs = generate_inputs_wannier90() + inputs["wannier90"]["wannier90"]["parameters"] = orm.Dict(w90_parameters) + + workchain = generate_workchain("wannier90_workflows.wannier90", inputs) + workchain.setup() + workchain.ctx.workchain_nscf = _generate_finished_nscf_calc( + generate_calc_job_node, fixture_localhost + ) + + return workchain + + +def test_prepare_wannier90_pp_inputs_fermi_from_nscf( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, +): # pylint: disable=redefined-outer-name + """Without a Fermi energy in the parameters, read it from the nscf run in the context.""" + workchain = _generate_workchain_without_scf_context( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, + w90_parameters={}, + ) + + inputs = workchain.prepare_wannier90_pp_inputs() + + assert inputs["wannier90"]["parameters"]["fermi_energy"] == pytest.approx( + NSCF_SCF_FERMI_ENERGY + ) + + +def test_prepare_wannier90_pp_inputs_fermi_from_parameters( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, +): # pylint: disable=redefined-outer-name + """A Fermi energy already in the parameters wins over the nscf run in the context.""" + given_fermi_energy = 1.23 + assert given_fermi_energy != NSCF_SCF_FERMI_ENERGY + + workchain = _generate_workchain_without_scf_context( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, + w90_parameters={"fermi_energy": given_fermi_energy}, + ) + + inputs = workchain.prepare_wannier90_pp_inputs() + + assert inputs["wannier90"]["parameters"]["fermi_energy"] == pytest.approx( + given_fermi_energy + ) + + +def test_prepare_wannier90_pp_inputs_rejects_none_fermi_in_parameters( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, +): # pylint: disable=redefined-outer-name + """A `fermi_energy` of None in the parameters is reported against the parameters.""" + workchain = _generate_workchain_without_scf_context( + generate_workchain, + generate_inputs_wannier90, + generate_calc_job_node, + fixture_localhost, + w90_parameters={"fermi_energy": None}, + ) + + with pytest.raises(ValueError, match="wannier90 parameters is None"): + workchain.prepare_wannier90_pp_inputs()