diff --git a/express/parsers/apps/nwchem/formats/txt.py b/express/parsers/apps/nwchem/formats/txt.py index c2886eaa..8262a0d3 100644 --- a/express/parsers/apps/nwchem/formats/txt.py +++ b/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. diff --git a/express/parsers/apps/nwchem/parser.py b/express/parsers/apps/nwchem/parser.py index 667af6f9..b34492e6 100644 --- a/express/parsers/apps/nwchem/parser.py +++ b/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): """ diff --git a/express/parsers/apps/nwchem/settings.py b/express/parsers/apps/nwchem/settings.py index 9494ca6d..b9f26d78 100644 --- a/express/parsers/apps/nwchem/settings.py +++ b/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 (?:(?PAlpha|Beta) )?Molecular Orbital Analysis") +VECTOR_REGEX = re.compile( + r"Vector\s+(?P\d+)\s+Occ=\s*(?P[\dDEe.+-]+)\s+E=\s*(?P[\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, diff --git a/express/parsers/mixins/electronic.py b/express/parsers/mixins/electronic.py index 1826c47a..65e5aba3 100644 --- a/express/parsers/mixins/electronic.py +++ b/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): """ diff --git a/express/parsers/utils.py b/express/parsers/utils.py index bf3312ea..63e1d895 100644 --- a/express/parsers/utils.py +++ b/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")) diff --git a/tests/fixtures/nwchem/references.py b/tests/fixtures/nwchem/references.py index 1db25356..abb10958 100644 --- a/tests/fixtures/nwchem/references.py +++ b/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 diff --git a/tests/integration/parsers/apps/nwchem/test_parser.py b/tests/integration/parsers/apps/nwchem/test_parser.py index 269a47a7..a0fbc66b 100644 --- a/tests/integration/parsers/apps/nwchem/test_parser.py +++ b/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) diff --git a/tests/manifest.yaml b/tests/manifest.yaml index f0bb42f2..aa9c6001 100644 --- a/tests/manifest.yaml +++ b/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 diff --git a/tests/unit/parsers/test_nwchem_txt_parser.py b/tests/unit/parsers/test_nwchem_txt_parser.py new file mode 100644 index 00000000..a3256d54 --- /dev/null +++ b/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"), [])