units as parameter in read_*

CHANGELOG.md

@@ -9,6 +9,7 @@ we hit release version 1.0.0.
 ## [0.15.0] - YYYY-MM-DD
 
 ### Added
+- `units` to `read_*` for some `Sile`s, #726
 - "Hz", "MHz", "GHz", "THz", and "invcm" as valid energy units, #725
 - added `read_gtensor` and `read_hyperfine_coupling` to `txtSileORCA`, #722
 - enabled `AtomsArgument` and `OrbitalsArgument` to accept `bool` for *all* or *none*

src/sisl/io/orca/stdout.py

@@ -5,7 +5,7 @@
 
 from sisl._internal import set_module
 from sisl.messages import deprecation
-from sisl.unit import units
+from sisl.unit import serialize_units_arg, unit_convert
 from sisl.utils import PropertyDict
 
 from .._multiple import SileBinder
@@ -260,12 +260,21 @@ def read_block(step_to):
 
     @SileBinder()
     @sile_fh_open()
-    def read_energy(self):
+    def read_energy(self, units="eV"):
         """Reads the energy blocks
 
+        Parameters
+        ----------
+        units : {str, dict, list, tuple}
+            selects units in the returned data
+
+        Note
+        ----
+        Energies written by ORCA have units of Ha.
+
         Returns
         -------
-        PropertyDict or list of PropertyDict : all energy data (in eV) from the "TOTAL SCF ENERGY" and "DFT DISPERSION CORRECTION" blocks
+        PropertyDict or list of PropertyDict : all energy data from the "TOTAL SCF ENERGY" and "DFT DISPERSION CORRECTION" blocks
         """
         f = self.step_to("TOTAL SCF ENERGY", allow_reread=False)[0]
         if not f:
@@ -274,28 +283,30 @@ def read_energy(self):
         self.readline()  # skip ---
         self.readline()  # skip blank line
 
-        Ha2eV = units("Ha", "eV")
+        units = serialize_units_arg(units)
+        Ha2unit = unit_convert("Ha", units["energy"])
+
         E = PropertyDict()
 
         line = self.readline()
         while "----" not in line:
             v = line.split()
             if "Total Energy" in line:
-                E["total"] = float(v[-4]) * Ha2eV
+                E["total"] = float(v[-4]) * Ha2unit
             elif "E(X)" in line:
-                E["exchange"] = float(v[-2]) * Ha2eV
+                E["exchange"] = float(v[-2]) * Ha2unit
             elif "E(C)" in line:
-                E["correlation"] = float(v[-2]) * Ha2eV
+                E["correlation"] = float(v[-2]) * Ha2unit
             elif "E(XC)" in line:
-                E["xc"] = float(v[-2]) * Ha2eV
+                E["xc"] = float(v[-2]) * Ha2unit
             elif "DFET-embed. en." in line:
-                E["embedding"] = float(v[-2]) * Ha2eV
+                E["embedding"] = float(v[-2]) * Ha2unit
             line = self.readline()
 
         if self.info.vdw_correction:
             self.step_to("DFT DISPERSION CORRECTION")
             v = self.step_to("Dispersion correction", allow_reread=False)[1].split()
-            E["vdw"] = float(v[-1]) * Ha2eV
+            E["vdw"] = float(v[-1]) * Ha2unit
 
         return E
 

src/sisl/io/orca/tests/test_txt.py

@@ -113,12 +113,12 @@ def test_hyperfine_coupling(sisl_files):
     # file without hyperfine_coupling tensors
     f = sisl_files(_dir, "nitric_oxide", "molecule_property.txt")
     out = txtSileORCA(f)
-    assert out.read_hyperfine_coupling() is None
+    assert out.read_hyperfine_coupling(units="MHz") is None
 
     # file with hyperfine_coupling tensors
     f = sisl_files(_dir, "phenalenyl", "molecule_property.txt")
     out = txtSileORCA(f)
-    A = out.read_hyperfine_coupling()
+    A = out.read_hyperfine_coupling(units="MHz")
     assert len(A) == 22
     assert A[0].iso == -23.380794
     assert A[1].ia == 1
@@ -136,3 +136,19 @@ def test_hyperfine_coupling(sisl_files):
     assert A[1].iso == 26.247902
     assert A[12].sa == "C"
     assert A[13].sa == "H"
+
+
+def test_hyperfine_coupling_units(sisl_files):
+    f = sisl_files(_dir, "phenalenyl", "molecule_property.txt")
+    out = txtSileORCA(f)
+    A = out.read_hyperfine_coupling()
+    B = out.read_hyperfine_coupling(units="eV")
+    C = out.read_hyperfine_coupling(units={"energy": "eV"})
+
+    assert A[0].iso == B[0].iso == C[0].iso
+
+    A = out.read_hyperfine_coupling(units=("MHz", "Ang"))
+    B = out.read_hyperfine_coupling(units={"energy": "MHz"})
+
+    assert A[0].iso == B[0].iso
+    assert A[0].iso != C[0].iso

src/sisl/io/orca/txt.py

@@ -6,7 +6,7 @@
 from sisl._core.geometry import Geometry
 from sisl._internal import set_module
 from sisl.messages import deprecation
-from sisl.unit import units
+from sisl.unit import serialize_units_arg, unit_convert
 from sisl.utils import PropertyDict
 
 from .._multiple import SileBinder
@@ -75,12 +75,21 @@ def read_electrons(self):
 
     @SileBinder()
     @sile_fh_open()
-    def read_energy(self):
+    def read_energy(self, units="eV"):
         """Reads the energy blocks
 
+        Parameters
+        ----------
+        units : {str, dict, list, tuple}
+            selects units in the returned data
+
+        Note
+        ----
+        Energies written by ORCA have units of Ha.
+
         Returns
         -------
-        PropertyDict or list of PropertyDict : all data (in eV) from the "DFT_Energy" and "VdW_Correction" blocks
+        PropertyDict or list of PropertyDict : all data from the "DFT_Energy" and "VdW_Correction" blocks
         """
         # read the DFT_Energy block
         f = self.step_to("$ DFT_Energy", allow_reread=False)[0]
@@ -90,13 +99,15 @@ def read_energy(self):
         self.readline()  # geom. index
         self.readline()  # prop. index
 
-        Ha2eV = units("Ha", "eV")
+        units = serialize_units_arg(units)
+        Ha2unit = unit_convert("Ha", units["energy"])
+
         E = PropertyDict()
 
         line = self.readline()
         while "----" not in line:
             v = line.split()
-            value = float(v[-1]) * Ha2eV
+            value = float(v[-1]) * Ha2unit
             if v[0] == "Exchange":
                 E["exchange"] = value
             elif v[0] == "Correlation":
@@ -115,7 +126,7 @@ def read_energy(self):
         line = self.readline()
         if self.info.vdw_correction:
             v = self.step_to("Van der Waals Correction:")[1].split()
-            E["vdw"] = float(v[-1]) * Ha2eV
+            E["vdw"] = float(v[-1]) * Ha2unit
 
         return E
 
@@ -190,9 +201,14 @@ def read_gtensor(self):
         return G
 
     @sile_fh_open()
-    def read_hyperfine_coupling(self):
+    def read_hyperfine_coupling(self, units="eV"):
         r"""Reads hyperfine couplings from the ``EPRNMR_ATensor`` block
 
+        Parameters
+        ----------
+        units : {str, dict, list, tuple}
+           selects units in the returned data
+
         For a nucleus :math:`k`, the hyperfine interaction is usually
         written in terms of the symmetric :math:`3\times 3` hyperfine
         tensor :math:`A_k` such that
@@ -209,8 +225,7 @@ def read_hyperfine_coupling(self):
 
         Note
         ----
-        The hyperfine tensors are given in units of MHz. This may change
-        in later versions.
+        Hyperfine tensors written by ORCA have units of MHz.
 
         Returns
         -------
@@ -220,6 +235,9 @@ def read_hyperfine_coupling(self):
         if not f:
             return None
 
+        units = serialize_units_arg(units)
+        MHz2unit = unit_convert("MHz", units["energy"])
+
         def read_A():
             A = PropertyDict()
             # Read the EPRNMR_ATensor block
@@ -238,7 +256,7 @@ def read_A():
             A["spin"] = float(v[-1])
 
             v = self.readline().split()
-            A["prefactor"] = float(v[-1])
+            A["prefactor"] = float(v[-1]) * MHz2unit
 
             self.readline()  # skip Raw line
             self.readline()  # skip header
@@ -247,7 +265,7 @@ def read_A():
             for i in range(3):
                 v = self.readline().split()
                 tensor[i] = v[1:4]
-            A["tensor"] = tensor  # raw A_total tensor
+            A["tensor"] = tensor * MHz2unit  # raw A_total tensor
 
             self.readline()  # skip eigenvector line
             self.readline()  # skip header
@@ -262,10 +280,12 @@ def read_A():
             self.readline()  # skip header
 
             v = self.readline().split()
-            A["eigenvalues"] = np.array(v[1:4], np.float64)  # eigenvalues of A_total
+            A["eigenvalues"] = (
+                np.array(v[1:4], np.float64) * MHz2unit
+            )  # eigenvalues of A_total
 
             v = self.readline().split()
-            A["iso"] = float(v[1])  # Fermi contact A_FC
+            A["iso"] = float(v[1]) * MHz2unit  # Fermi contact A_FC
 
             return A
 

src/sisl/io/vasp/doscar.py

@@ -5,6 +5,7 @@
 
 from sisl._array import arrayf
 from sisl._internal import set_module
+from sisl.unit import serialize_units_arg, unit_convert
 
 from ..sile import add_sile, sile_fh_open
 from .sile import SileVASP
@@ -33,15 +34,20 @@ def read_fermi_level(self):
         return float(line[3])
 
     @sile_fh_open()
-    def read_data(self):
+    def read_data(self, units="eV"):
         r"""Read DOS, as calculated and written by VASP
 
+        Parameters
+        ----------
+        units : {str, dict, list, tuple}
+            selects units in the returned data
+
         Returns
         -------
         E : numpy.ndarray
-            energy points (in eV)
+            energy points
         DOS : numpy.ndarray
-            DOS points (in 1/eV)
+            DOS points (units of 1/energy)
         """
         # read first line
         self.readline()  # NIONS, NIONS, JOBPAR_, WDES%INCDIJ
@@ -66,7 +72,11 @@ def read_data(self):
             line = arrayf(self.readline().split())
             E[ie] = line[0]
             DOS[:, ie] = line[1 : ns + 1]
-        return E - Ef, DOS
+
+        units = serialize_units_arg(units)
+        eV2unit = unit_convert("eV", units["energy"])
+
+        return (E - Ef) * eV2unit, DOS / eV2unit
 
 
 add_sile("DOSCAR", doscarSileVASP, gzip=True)

src/sisl/io/vasp/eigenval.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from sisl._internal import set_module
+from sisl.unit import serialize_units_arg, unit_convert
 
 from ..sile import add_sile, sile_fh_open
 
@@ -18,13 +19,15 @@ class eigenvalSileVASP(SileVASP):
     """Kohn-Sham eigenvalues"""
 
     @sile_fh_open()
-    def read_data(self, k=False):
+    def read_data(self, k=False, units="eV"):
         r"""Read eigenvalues, as calculated and written by VASP
 
         Parameters
         ----------
         k : bool, optional
            also return k points and weights
+        units : {str, dict, list, tuple}
+           selects units in the returned data
 
         Returns
         -------
@@ -58,6 +61,11 @@ def read_data(self, k=False):
                 # We currently neglect the populations
                 E = map(float, self.readline().split()[1 : ns + 1])
                 eigs[:, ik, ib] = list(E)
+
+        units = serialize_units_arg(units)
+        eV2unit = unit_convert("eV", units["energy"])
+        eigs *= eV2unit
+
         if k:
             return eigs, kk, w
         return eigs

src/sisl/io/vasp/locpot.py

@@ -7,6 +7,7 @@
 
 from sisl import Grid
 from sisl._internal import set_module
+from sisl.unit import serialize_units_arg, unit_convert
 
 from .._help import grid_reduce_indices
 from ..sile import add_sile, sile_fh_open
@@ -24,8 +25,8 @@ class locpotSileVASP(carSileVASP):
     """
 
     @sile_fh_open(True)
-    def read_grid(self, index=0, dtype=np.float64, **kwargs):
-        """Reads the potential (in eV) from the file and returns with a grid (plus geometry)
+    def read_grid(self, index=0, dtype=np.float64, units="eV", **kwargs):
+        """Reads the potential from the file and returns with a grid (plus geometry)
 
         Parameters
         ----------
@@ -37,6 +38,8 @@ def read_grid(self, index=0, dtype=np.float64, **kwargs):
            contributions for each corresponding index.
         dtype : numpy.dtype, optional
            grid stored dtype
+        units : {str, dict, list, tuple}
+           selects units in the returned data
         spin : optional
            same as `index` argument. `spin` argument has precedence.
 
@@ -48,6 +51,9 @@ def read_grid(self, index=0, dtype=np.float64, **kwargs):
         geom = self.read_geometry()
         V = geom.lattice.volume
 
+        units = serialize_units_arg(units)
+        eV2unit = unit_convert("eV", units["energy"])
+
         # Now we are past the cell and geometry
         # We can now read the size of LOCPOT
         self.readline()
@@ -96,7 +102,7 @@ def read_grid(self, index=0, dtype=np.float64, **kwargs):
         # Since we populate the grid data afterwards there
         # is no need to create a bigger grid than necessary.
         grid = Grid([1, 1, 1], dtype=dtype, geometry=geom)
-        grid.grid = val
+        grid.grid = val * eV2unit
 
         return grid