💥 Units: switch to QE-native units for outputs

docs/getting_started.md

@@ -37,7 +37,7 @@ pw_out.get_output('fermi_energy')
 
 ### Converting to other units
 
-By default, `qe-tools` returns energies in eV. You can obtain a [`pint`](https://pint.readthedocs.io/en/stable/) quantity with unit attached using the `to` input (requires the `pint` extra: `pip install qe-tools[pint]`):
+By default, `pw.x` outputs are reported in QE-native units (Ry, bohr, kbar, ...). You can obtain a [`pint`](https://pint.readthedocs.io/en/stable/) quantity with unit attached using the `to` input (requires the `pint` extra: `pip install qe-tools[pint]`):
 
 ```python
 fermi_energy = pw_out.get_output('fermi_energy', to='pint')
@@ -47,10 +47,10 @@ fermi_energy
 and can then convert the value to any unit you prefer:
 
 ```python
-fermi_energy.to('Ry')
+fermi_energy.to('eV')
 ```
 
-See the [units section](units.md) for more information on the list of units we return quantities in by default.
+See the [units section](units.md) for the default units of each output class.
 
 ### Tab completion
 

docs/units.md

@@ -1,20 +1,29 @@
 # Units
 
-Quantum ESPRESSO uses different units for the same quantity across input, stdout, and XML.
-Instead, `qe-tools` exposes outputs on e.g. the `PwOutput` class in a single opinionated set of units, independent of what Quantum ESPRESSO writes to the XML or stdout.
-Below you can find the units we chose for each quantity:
-
-| Quantity              | Unit              |
-| --------------------- | ----------------- |
-| Energy                | `eV`              |
-| Force                 | `eV/Å`            |
-| Stress / pressure     | `GPa`             |
-| Length                | `Å`               |
-| Volume                | `ų`              |
-| Reciprocal length     | `1/Å`             |
-| Density of states     | `1/eV`            |
-| Magnetic moment       | `μ_B`             |
-| Time                  | `s`               |
+Currently, `qe-tools` aims to report each output in the unit Quantum ESPRESSO uses, rather than imposing a single unit choice across the board.
+This keeps the values you get out aligned with what most Quantum ESPRESSO users would expect.
+
+The unit of each output is declared on each field (via a `Unit(...)` marker) and stated explicitly in the field's docstring.
+The defaults per output class are:
+
+| Output class       | Energy | Length  | Force     | Stress | Other            |
+| ------------------ | ------ | ------- | --------- | ------ | ---------------- |
+| `PwOutput`         | `Ry`   | `bohr`  | `Ry/bohr` | `kbar` | volume `bohr³`, k-points `1/bohr`, magnetization `μ_B`, time `s` |
+| `DosOutput`        | `eV`   | —       | —         | —      | DOS `1/eV`       |
+| `BandsOutput`      | `eV`   | —       | —         | —      | k-points in the units `bands.x` writes (typically `2π/alat`) |
+| `ProjwfcOutput`    | `eV`   | —       | —         | —      | pDOS `1/eV`      |
+
+A couple of fields in `PwOutput` come from sources QE prints in eV (e.g. `highest_occupied_level`, `lowest_unoccupied_level` from the stdout line `... (ev): ...`); those are exposed in eV, as written.
+The docstring of every output field states its unit, so the source of truth is always the field itself.
+
+Some outputs - like the total energy - are written in Ry in the stdout but in Ha in the XML.
+Here we again choose the units most users will be familiar with: Ry.
+
+!!! info "Switching to a single set of units"
+
+    Originally, we had selected a single set of units (eV, eV/Å, ...) and wanted to enforce each output (and future input) was specified in these units.
+    However, this turned out to cause a lot of friction with the QE ecosystem.
+    We're looking into making the unit set configurable, so this consistent behavior can be recovered in the future.
 
 ## Converting to different units
 
@@ -29,9 +38,9 @@ You can then request a `pint` quantity from any unit-annotated output:
 fermi_energy = pw_out.get_output('fermi_energy', to='pint')
 ```
 
-This will return a `pint` quantity with unit.
+This will return a `pint` quantity with unit attached.
 You can then convert the value to any unit you prefer:
 
 ```python
-fermi_energy.to('Ry')
+fermi_energy.to('eV')
 ```

src/qe_tools/converters/aiida.py

@@ -12,6 +12,8 @@ class AiiDAConverter(BaseConverter):
     def get_conversion_mapping(cls) -> dict[str, typing.Any]:
         import numpy as np
 
+        from qe_tools import CONSTANTS
+
         try:
             from aiida import orm
         except ImportError:
@@ -45,8 +47,14 @@ def convert_dos(energy: list, dos: list | dict):
                 convert_structure_data,
                 {
                     "symbols": "atomic_species",
-                    "cell": ("cell", lambda cell: np.array(cell)),
-                    "positions": ("positions", lambda positions: np.array(positions)),
+                    "cell": (
+                        "cell",
+                        lambda cell: np.array(cell) * CONSTANTS.bohr_to_ang,
+                    ),
+                    "positions": (
+                        "positions",
+                        lambda positions: np.array(positions) * CONSTANTS.bohr_to_ang,
+                    ),
                 },
             ),
             "full_dos": (

src/qe_tools/converters/ase.py

@@ -10,6 +10,8 @@ class ASEConverter(BaseConverter):
     def get_conversion_mapping(cls) -> dict[str, typing.Any]:
         import numpy as np
 
+        from qe_tools import CONSTANTS
+
         try:
             from ase import Atoms
         except ImportError:
@@ -24,8 +26,14 @@ def get_conversion_mapping(cls) -> dict[str, typing.Any]:
                 Atoms,
                 {
                     "symbols": "symbols",
-                    "cell": ("cell", lambda cell: np.array(cell)),
-                    "positions": ("positions", lambda positions: np.array(positions)),
+                    "cell": (
+                        "cell",
+                        lambda cell: np.array(cell) * CONSTANTS.bohr_to_ang,
+                    ),
+                    "positions": (
+                        "positions",
+                        lambda positions: np.array(positions) * CONSTANTS.bohr_to_ang,
+                    ),
                 },
             ),
         }

src/qe_tools/converters/pymatgen.py

@@ -10,6 +10,8 @@ class PymatgenConverter(BaseConverter):
     def get_conversion_mapping(cls) -> dict[str, typing.Any]:
         import numpy as np
 
+        from qe_tools import CONSTANTS
+
         try:
             from pymatgen.core.structure import Structure
         except ImportError:
@@ -24,8 +26,14 @@ def get_conversion_mapping(cls) -> dict[str, typing.Any]:
                 Structure,
                 {
                     "species": "symbols",
-                    "lattice": ("cell", lambda cell: np.array(cell)),
-                    "coords": ("positions", lambda positions: np.array(positions)),
+                    "lattice": (
+                        "cell",
+                        lambda cell: np.array(cell) * CONSTANTS.bohr_to_ang,
+                    ),
+                    "coords": (
+                        "positions",
+                        lambda positions: np.array(positions) * CONSTANTS.bohr_to_ang,
+                    ),
                 },
             ),
         }

src/qe_tools/outputs/parsers/pw.py

@@ -125,17 +125,15 @@ def parse(content: str) -> dict:
             ]
 
         # k-points: take the last block, since vc-relax reprints after relaxation.
-        # The stdout values are in 2pi/alat; convert to 1/Å using `alat` (also from
+        # The stdout values are in 2pi/alat; convert to 1/bohr using `alat` (also from
         # stdout, in bohr) so downstream Specs match the XML-derived units.
         kpoint_blocks = list(_KPOINTS_BLOCK_RE.finditer(content))
         alat_match = _ALAT_RE.search(content)
         if kpoint_blocks and alat_match:
             import math
 
-            from qe_tools import CONSTANTS
-
-            alat_ang = float(alat_match.group(1)) * CONSTANTS.bohr_to_ang
-            scale = 2 * math.pi / alat_ang
+            alat_bohr = float(alat_match.group(1))
+            scale = 2 * math.pi / alat_bohr
             block = kpoint_blocks[-1]
             cartesian = []
             weights = []