update: read homo/lumo from converged values

express/parsers/apps/nwchem/formats/txt.py

@@ -1,6 +1,7 @@
 from express.parsers.settings import Constant  # noqa: F401
 from express.parsers.apps.nwchem import settings
 from express.parsers.formats.txt import BaseTXTParser
+from express.parsers.utils import _fortran_float
 
 
 class NwchemTXTParser(BaseTXTParser):
@@ -11,6 +12,34 @@ class NwchemTXTParser(BaseTXTParser):
     def __init__(self, work_dir):
         super(NwchemTXTParser, self).__init__(work_dir)
 
+    def eigenvalues_at_vectors(self, text):
+        """
+        Extracts eigenvalues at molecular orbitals (vectors). Geometry optimizations print one
+        orbital analysis section per step; the last one is read, for the final geometry. A
+        spin-polarized run prints an alpha and a beta section per step, both of which are read.
+
+        Units:
+            energy: Hartree
+
+        Args:
+            text (str): text to extract data from.
+
+        Returns:
+            list[dict]
+        """
+        blocks = list(settings.ORBITAL_ANALYSIS_BLOCK_REGEX.finditer(text))
+        ends = [block.start() for block in blocks[1:]] + [len(text)]
+        last_block_per_spin = {block.group("spin"): text[block.end() : end] for block, end in zip(blocks, ends)}
+        return [
+            {
+                "vector": int(orbital.group("vector")),
+                "occupation": _fortran_float(orbital.group("occupation")),
+                "energy": _fortran_float(orbital.group("energy")),
+            }
+            for block in last_block_per_spin.values()
+            for orbital in settings.VECTOR_REGEX.finditer(block)
+        ]
+
     def total_energy(self, text):
         """
         Extracts total energy.
@@ -40,30 +69,6 @@ def total_energy_contributions(self, text):
                 energy_contributions.update({contribution: {"name": contribution, "value": value}})
         return energy_contributions
 
-    def homo_energy(self, text):
-        """
-        Extracts HOMO energy.
-
-        Args:
-            text (str): text to extract data from.
-
-        Returns:
-            float | None
-        """
-        return self._general_output_parser(text, **settings.REGEX["homo_energy"])
-
-    def lumo_energy(self, text):
-        """
-        Extracts LUMO energy.
-
-        Args:
-            text (str): text to extract data from.
-
-        Returns:
-            float | None
-        """
-        return self._general_output_parser(text, **settings.REGEX["lumo_energy"])
-
     def zero_point_energy(self, text):
         """
         Extracts zero point energy.

express/parsers/apps/nwchem/parser.py

@@ -49,25 +49,30 @@ def total_energy_contributions(self):
                     value1[key2] = value2 * Constant.HARTREE
         return energy_contributions
 
-    def homo_energy(self):
+    def eigenvalues_at_vectors(self):
         """
-        Returns HOMO energy.
+        Returns eigenvalues at molecular orbitals (vectors).
 
         Reference:
+            func: express.parsers.mixins.electronic.ElectronicDataMixin.eigenvalues_at_vectors
             NWChem orbital energies are defaulted to hartrees and are converted to eV in this method.
         """
-        homo_energy = self.txt_parser.homo_energy(self._get_file_content(self.stdout_file))
-        return None if homo_energy is None else Constant.HARTREE * homo_energy
+        orbitals = self.txt_parser.eigenvalues_at_vectors(self._get_file_content(self.stdout_file))
+        return [dict(orbital, energy=Constant.HARTREE * orbital["energy"]) for orbital in orbitals]
 
-    def lumo_energy(self):
+    def homo_energy(self):
+        """
+        Returns HOMO energy, the highest energy among the occupied molecular orbitals.
         """
-        Returns LUMO energy.
+        energies = [orbital["energy"] for orbital in self.eigenvalues_at_vectors() if orbital["occupation"] > 0]
+        return max(energies) if energies else None
 
-        Reference:
-            NWChem orbital energies are defaulted to hartrees and are converted to eV in this method.
+    def lumo_energy(self):
+        """
+        Returns LUMO energy, the lowest energy among the unoccupied molecular orbitals.
         """
-        lumo_energy = self.txt_parser.lumo_energy(self._get_file_content(self.stdout_file))
-        return None if lumo_energy is None else Constant.HARTREE * lumo_energy
+        energies = [orbital["energy"] for orbital in self.eigenvalues_at_vectors() if orbital["occupation"] == 0]
+        return min(energies) if energies else None
 
     def zero_point_energy(self):
         """

express/parsers/apps/nwchem/settings.py

@@ -1,13 +1,19 @@
+import re
+
 from express.parsers.settings import GENERAL_REGEX
 
 COMMON_REGEX = r"{}\s+[=:<>]\s*([-+]?\d*\.?\d*([Ee][+-]?\d+)?)"
 DOUBLE_REGEX = GENERAL_REGEX["double_number"]
 NWCHEM_OUTPUT_FILE_REGEX = "Northwest Computational Chemistry Package"
 
+# Closed-shell runs print a single unlabeled section, spin-polarized (ODFT) ones an Alpha and a Beta.
+ORBITAL_ANALYSIS_BLOCK_REGEX = re.compile(r"DFT Final (?:(?P<spin>Alpha|Beta) )?Molecular Orbital Analysis")
+VECTOR_REGEX = re.compile(
+    r"Vector\s+(?P<vector>\d+)\s+Occ=\s*(?P<occupation>[\dDEe.+-]+)\s+E=\s*(?P<energy>[\dDEe.+-]+)"
+)
+
 REGEX = {
     "total_energy": {"regex": COMMON_REGEX.format("Total DFT energy"), "occurrences": -1, "output_type": "float"},
-    "homo_energy": {"regex": COMMON_REGEX.format("HOMO"), "occurrences": -1, "output_type": "float"},
-    "lumo_energy": {"regex": COMMON_REGEX.format("LUMO"), "occurrences": -1, "output_type": "float"},
     "zero_point_energy": {
         "regex": COMMON_REGEX.format("Zero-Point correction to Energy"),
         "occurrences": -1,

express/parsers/mixins/electronic.py

@@ -85,6 +85,28 @@ def eigenvalues_at_kpoints(self):
         """
         pass
 
+    @abstractmethod
+    def eigenvalues_at_vectors(self):
+        """
+        Returns eigenvalues at molecular orbitals (vectors), the molecular analogue of
+        `eigenvalues_at_kpoints` for systems without reciprocal space.
+
+        Units:
+            energy: eV
+
+        Returns:
+             list[dict]
+
+        Example:
+            [
+                {'vector': 1, 'occupation': 2.0, 'energy': -520.7098},
+                {'vector': 2, 'occupation': 2.0, 'energy': -27.2698},
+                {'vector': 3, 'occupation': 0.0, 'energy': 1.7931},
+                ...
+            ]
+        """
+        pass
+
     @abstractmethod
     def dos(self):
         """

express/parsers/utils.py

@@ -80,3 +80,18 @@ def lattice_basis_to_poscar(lattice: dict, basis: dict, basis_units: str = "cart
             "\n".join([" ".join(["{0:14.9f}".format(v) for v in x["value"]]) for x in basis["coordinates"]]),
         ]
     )
+
+
+
+def _fortran_float(value):
+    """
+    Converts a Fortran-formatted float string (e.g. "-1.234D+01", double-precision
+    "D" exponent notation) to a Python float.
+
+    Args:
+        value (str): Fortran-formatted number, e.g. "-1.234D+01".
+
+    Returns:
+        float
+    """
+    return float(value.replace("D", "E").replace("d", "e"))

tests/fixtures/nwchem/references.py

@@ -3,9 +3,31 @@
 All nwchem output values are in hartrees. ExPrESS converts units to eV.
 All reference energies are in eV.
 """
+
 TOTAL_ENERGY = -2079.18666382721904
-HOMO_ENERGY = -12.800485418916242
-LUMO_ENERGY = 3.1242763921882197
+
+# Only the count and the edges of the spectrum are pinned; HOMO/LUMO below cover the frontier.
+EIGENVALUES_AT_VECTORS_COUNT = 19
+EIGENVALUES_AT_VECTORS_FIRST = {"vector": 1, "occupation": 2.0, "energy": -520.7096480731956}
+EIGENVALUES_AT_VECTORS_LAST = {"vector": 19, "occupation": 0.0, "energy": 96.85556141135692}
+
+# The *_INITIAL_GUESS values come from the initial-guess section preceding the SCF cycle; the tests
+# assert the parser does not return them.
+HOMO_ENERGY = -7.938587261191046
+LUMO_ENERGY = 1.793148251055798
+HOMO_ENERGY_INITIAL_GUESS = -12.800485418916242
+LUMO_ENERGY_INITIAL_GUESS = 3.1242763921882197
+
+# test-002/nwchem-frequency.log holds one orbital analysis section per optimization step.
+EIGENVALUES_AT_VECTORS_MULTISTEP_COUNT = 18
+EIGENVALUES_AT_VECTORS_MULTISTEP_FIRST = {"vector": 1, "occupation": 2.0, "energy": -520.5444749015667}
+EIGENVALUES_AT_VECTORS_MULTISTEP_LAST = {"vector": 18, "occupation": 0.0, "energy": 70.15805296279406}
+
+HOMO_ENERGY_MULTISTEP = -7.8966272890902385
+LUMO_ENERGY_MULTISTEP = 1.7759215328081597
+HOMO_ENERGY_MULTISTEP_INITIAL_GUESS = -10.187320671325603
+LUMO_ENERGY_MULTISTEP_INITIAL_GUESS = -3.5424611206222094
+
 ZERO_POINT_ENERGY = 0.5748347036575007
 THERMAL_CORRECTION_TO_ENERGY = 15.033
 THERMAL_CORRECTION_TO_ENTHALPY = 15.626

tests/integration/parsers/apps/nwchem/test_parser.py

@@ -15,11 +15,37 @@ def tearDown(self):
     def test_nwchem_total_energy(self):
         self.assertAlmostEqual(self.parser.total_energy(), TOTAL_ENERGY, places=2)
 
+    def test_nwchem_eigenvalues_at_vectors(self):
+        eigenvalues = self.parser.eigenvalues_at_vectors()
+        self.assertEqual(len(eigenvalues), EIGENVALUES_AT_VECTORS_COUNT)
+        self.assertDeepAlmostEqual(eigenvalues[0], EIGENVALUES_AT_VECTORS_FIRST, places=2)
+        self.assertDeepAlmostEqual(eigenvalues[-1], EIGENVALUES_AT_VECTORS_LAST, places=2)
+
+    def test_nwchem_eigenvalues_at_vectors_multistep(self):
+        eigenvalues = self.parser.eigenvalues_at_vectors()
+        self.assertEqual(len(eigenvalues), EIGENVALUES_AT_VECTORS_MULTISTEP_COUNT)
+        self.assertDeepAlmostEqual(eigenvalues[0], EIGENVALUES_AT_VECTORS_MULTISTEP_FIRST, places=2)
+        self.assertDeepAlmostEqual(eigenvalues[-1], EIGENVALUES_AT_VECTORS_MULTISTEP_LAST, places=2)
+
     def test_nwchem_homo_energy(self):
-        self.assertAlmostEqual(self.parser.homo_energy(), HOMO_ENERGY, places=2)
+        homo_energy = self.parser.homo_energy()
+        self.assertAlmostEqual(homo_energy, HOMO_ENERGY, places=2)
+        self.assertNotAlmostEqual(homo_energy, HOMO_ENERGY_INITIAL_GUESS, places=2)
 
     def test_nwchem_lumo_energy(self):
-        self.assertAlmostEqual(self.parser.lumo_energy(), LUMO_ENERGY, places=2)
+        lumo_energy = self.parser.lumo_energy()
+        self.assertAlmostEqual(lumo_energy, LUMO_ENERGY, places=2)
+        self.assertNotAlmostEqual(lumo_energy, LUMO_ENERGY_INITIAL_GUESS, places=2)
+
+    def test_nwchem_homo_energy_multistep(self):
+        homo_energy = self.parser.homo_energy()
+        self.assertAlmostEqual(homo_energy, HOMO_ENERGY_MULTISTEP, places=2)
+        self.assertNotAlmostEqual(homo_energy, HOMO_ENERGY_MULTISTEP_INITIAL_GUESS, places=2)
+
+    def test_nwchem_lumo_energy_multistep(self):
+        lumo_energy = self.parser.lumo_energy()
+        self.assertAlmostEqual(lumo_energy, LUMO_ENERGY_MULTISTEP, places=2)
+        self.assertNotAlmostEqual(lumo_energy, LUMO_ENERGY_MULTISTEP_INITIAL_GUESS, places=2)
 
     def test_nwchem_total_energy_contributions(self):
         self.assertDeepAlmostEqual(self.parser.total_energy_contributions(), TOTAL_ENERGY_CONTRIBUTION, places=2)

tests/manifest.yaml

@@ -2,6 +2,14 @@ test_nwchem_total_energy:
   workDir: fixtures/nwchem/test-001
   stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log
 
+test_nwchem_eigenvalues_at_vectors:
+  workDir: fixtures/nwchem/test-001
+  stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log
+
+test_nwchem_eigenvalues_at_vectors_multistep:
+  workDir: fixtures/nwchem/test-002
+  stdoutFile: fixtures/nwchem/test-002/nwchem-frequency.log
+
 test_nwchem_homo_energy:
   workDir: fixtures/nwchem/test-001
   stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log
@@ -10,6 +18,14 @@ test_nwchem_lumo_energy:
   workDir: fixtures/nwchem/test-001
   stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log
 
+test_nwchem_homo_energy_multistep:
+  workDir: fixtures/nwchem/test-002
+  stdoutFile: fixtures/nwchem/test-002/nwchem-frequency.log
+
+test_nwchem_lumo_energy_multistep:
+  workDir: fixtures/nwchem/test-002
+  stdoutFile: fixtures/nwchem/test-002/nwchem-frequency.log
+
 test_nwchem_total_energy_contributions:
   workDir: fixtures/nwchem/test-001
   stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log

tests/unit/parsers/test_nwchem_txt_parser.py

@@ -0,0 +1,60 @@
+import unittest
+
+from express.parsers.apps.nwchem.formats.txt import NwchemTXTParser
+
+
+def orbital_analysis_block(spin, energies):
+    header = f"DFT Final {spin + ' ' if spin else ''}Molecular Orbital Analysis"
+    vectors = "\n".join(
+        f" Vector {index + 1:4d}  Occ={occupation}  E={energy}" for index, (occupation, energy) in enumerate(energies)
+    )
+    return f"{header}\n{'-' * len(header)}\n\n{vectors}\n\n center of mass\n"
+
+
+class NwchemTXTParserTest(unittest.TestCase):
+    """
+    Covers the selection of orbital analysis sections, which the integration fixtures cannot
+    exercise: both of them are closed shell and single channel.
+    """
+
+    def setUp(self):
+        self.parser = NwchemTXTParser(work_dir=".")
+
+    def test_reads_last_section_of_closed_shell_output(self):
+        text = orbital_analysis_block(None, [("2.000000D+00", "-1.0D+00")]) + orbital_analysis_block(
+            None, [("2.000000D+00", "-2.0D+00"), ("0.000000D+00", "3.0D-01")]
+        )
+        self.assertEqual(
+            self.parser.eigenvalues_at_vectors(text),
+            [
+                {"vector": 1, "occupation": 2.0, "energy": -2.0},
+                {"vector": 2, "occupation": 0.0, "energy": 0.3},
+            ],
+        )
+
+    def test_reads_both_channels_of_spin_polarized_output(self):
+        text = orbital_analysis_block("Alpha", [("1.000000D+00", "-9.0D+00")]) + orbital_analysis_block(
+            "Beta", [("1.000000D+00", "-8.0D+00"), ("0.000000D+00", "5.0D-01")]
+        )
+        self.assertEqual(
+            self.parser.eigenvalues_at_vectors(text),
+            [
+                {"vector": 1, "occupation": 1.0, "energy": -9.0},
+                {"vector": 1, "occupation": 1.0, "energy": -8.0},
+                {"vector": 2, "occupation": 0.0, "energy": 0.5},
+            ],
+        )
+
+    def test_reads_last_step_of_each_channel(self):
+        text = "".join(
+            [
+                orbital_analysis_block("Alpha", [("1.000000D+00", "-9.0D+00")]),
+                orbital_analysis_block("Beta", [("1.000000D+00", "-8.0D+00")]),
+                orbital_analysis_block("Alpha", [("1.000000D+00", "-7.0D+00")]),
+                orbital_analysis_block("Beta", [("1.000000D+00", "-6.0D+00")]),
+            ]
+        )
+        self.assertEqual([orbital["energy"] for orbital in self.parser.eigenvalues_at_vectors(text)], [-7.0, -6.0])
+
+    def test_returns_nothing_without_an_orbital_analysis_section(self):
+        self.assertEqual(self.parser.eigenvalues_at_vectors(" Total DFT energy = -76.4\n"), [])