diff --git a/express/parsers/__init__.py b/express/parsers/__init__.py index d2e233c3..03e85954 100644 --- a/express/parsers/__init__.py +++ b/express/parsers/__init__.py @@ -8,6 +8,8 @@ class BaseParser(RoundNumericValuesMixin): Base Parser class. """ + is_non_periodic = False + def __init__(self, *args, **kwargs): self.args = args self.kwargs = kwargs diff --git a/express/parsers/apps/nwchem/formats/txt.py b/express/parsers/apps/nwchem/formats/txt.py index 8262a0d3..65a7ad47 100644 --- a/express/parsers/apps/nwchem/formats/txt.py +++ b/express/parsers/apps/nwchem/formats/txt.py @@ -1,4 +1,4 @@ -from express.parsers.settings import Constant # noqa: F401 +from express.parsers.settings import Constant from express.parsers.apps.nwchem import settings from express.parsers.formats.txt import BaseTXTParser from express.parsers.utils import _fortran_float @@ -12,6 +12,53 @@ class NwchemTXTParser(BaseTXTParser): def __init__(self, work_dir): super(NwchemTXTParser, self).__init__(work_dir) + def geometry_blocks(self, text): + """ + Extracts every "Output coordinates" block as elements and angstrom coordinates -- the header + carries the a.u. factor, so `units au` runs read correctly too -- or nothing at all if the + blocks disagree on their atoms: a log truncated mid-table parses as a shorter molecule, and + rupy would publish the fragment as `final_structure`. + """ + blocks = [] + for block in settings.GEOMETRY_BLOCK_REGEX.finditer(text): + to_angstrom = float(block.group("scale")) * Constant.BOHR + rows = list(settings.GEOMETRY_ROW_REGEX.finditer(block.group("rows"))) + elements = [row.group("element") for row in rows] + blocks.append((elements, [[float(row.group(axis)) * to_angstrom for axis in "xyz"] for row in rows])) + if len({tuple(elements) for elements, _ in blocks}) > 1: + return [] + return blocks + + def basis(self, text, index): + """ + Returns the geometry block at the given index. A single-point run prints exactly one block, + so index 0 and index -1 then coincide. + + Args: + text (str): text to extract data from. + index (int): position of the block among those printed. + + Returns: + dict | None + + Example: + { + 'units': 'angstrom', + 'elements': [{'id': 0, 'value': 'O'}, {'id': 1, 'value': 'H'}], + 'coordinates': [{'id': 0, 'value': [0.0, 0.0, 0.11]}, {'id': 1, 'value': [0.0, 0.75, -0.46]}] + } + """ + blocks = self.geometry_blocks(text) + if not blocks: + return None + + elements, coordinates = blocks[index] + return { + "units": "angstrom", + "elements": [{"id": idx, "value": value} for idx, value in enumerate(elements)], + "coordinates": [{"id": idx, "value": coordinate} for idx, coordinate in enumerate(coordinates)], + } + def eigenvalues_at_vectors(self, text): """ Extracts eigenvalues at molecular orbitals (vectors). Geometry optimizations print one diff --git a/express/parsers/apps/nwchem/parser.py b/express/parsers/apps/nwchem/parser.py index b34492e6..6f3e5511 100644 --- a/express/parsers/apps/nwchem/parser.py +++ b/express/parsers/apps/nwchem/parser.py @@ -11,8 +11,13 @@ class NwchemParser(BaseParser, IonicDataMixin, ElectronicDataMixin, ReciprocalDataMixin): """ Nwchem parser class. + + NWChem prints no cell for a molecule, so the lattice methods are left unimplemented and return + None: `express.properties.utils.box_molecule` derives the box, rather than a parser inventing it. """ + is_non_periodic = True + def __init__(self, *args, **kwargs): super(NwchemParser, self).__init__(*args, **kwargs) self.work_dir = self.kwargs["work_dir"] @@ -49,6 +54,24 @@ def total_energy_contributions(self): value1[key2] = value2 * Constant.HARTREE return energy_contributions + def initial_basis(self): + """ + Returns initial basis. + + Reference: + func: express.parsers.mixins.ionic.IonicDataMixin.initial_basis + """ + return self.txt_parser.basis(self._get_file_content(self.stdout_file), 0) + + def final_basis(self): + """ + Returns final basis. + + Reference: + func: express.parsers.mixins.ionic.IonicDataMixin.final_basis + """ + return self.txt_parser.basis(self._get_file_content(self.stdout_file), -1) + def eigenvalues_at_vectors(self): """ Returns eigenvalues at molecular orbitals (vectors). diff --git a/express/parsers/apps/nwchem/settings.py b/express/parsers/apps/nwchem/settings.py index b9f26d78..0d130700 100644 --- a/express/parsers/apps/nwchem/settings.py +++ b/express/parsers/apps/nwchem/settings.py @@ -6,6 +6,39 @@ DOUBLE_REGEX = GENERAL_REGEX["double_number"] NWCHEM_OUTPUT_FILE_REGEX = "Northwest Computational Chemistry Package" +GEOMETRY_HEADER_REGEX = r"Output coordinates in (?P\S+) \(scale by\s+(?P[\d.]+) to convert to a\.u\.\)" +# Requiring the tag to start with a symbol keeps a numeric-looking row of another table from matching. +ELEMENT_REGEX = r"[A-Za-z]+" +GEOMETRY_TAG_REGEX = r"{}\S*".format(ELEMENT_REGEX) +GEOMETRY_RULE_REGEX = r"^[ \t]*-{4,}[- \t]*$\n" +GEOMETRY_ROW_TEMPLATE = r"^[ \t]*\d+[ \t]+{tag}[ \t]+{double}[ \t]+{x}[ \t]+{y}[ \t]+{z}[ \t]*$" +# Never across the next block, so a header is always paired with its own table. +UNTIL_NEXT_GEOMETRY_REGEX = r"(?:(?!Output coordinates in)[\s\S])*?" + +# The trailing `(?:row)+` bounds the table: it stops at the blank line after the last atom. +GEOMETRY_BLOCK_REGEX = re.compile( + GEOMETRY_HEADER_REGEX + + UNTIL_NEXT_GEOMETRY_REGEX + + GEOMETRY_RULE_REGEX + + r"(?P(?:{})+)".format( + GEOMETRY_ROW_TEMPLATE.format( + tag=GEOMETRY_TAG_REGEX, double=DOUBLE_REGEX, x=DOUBLE_REGEX, y=DOUBLE_REGEX, z=DOUBLE_REGEX + ) + + r"\n" + ), + re.MULTILINE, +) +GEOMETRY_ROW_REGEX = re.compile( + GEOMETRY_ROW_TEMPLATE.format( + tag=r"(?P{})\S*".format(ELEMENT_REGEX), + double=DOUBLE_REGEX, + x=r"(?P{})".format(DOUBLE_REGEX), + y=r"(?P{})".format(DOUBLE_REGEX), + z=r"(?P{})".format(DOUBLE_REGEX), + ), + re.MULTILINE, +) + # 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( diff --git a/express/properties/material.py b/express/properties/material.py index c9b6c151..23f12dfc 100644 --- a/express/properties/material.py +++ b/express/properties/material.py @@ -13,6 +13,7 @@ from express.properties.scalar.volume import Volume from express.properties.structural.inchi import Inchi from express.properties.structural.inchi_key import InchiKey +from express.properties.utils import box_molecule class Material(BaseProperty): @@ -22,7 +23,7 @@ class Material(BaseProperty): def __init__(self, name, parser, *args, **kwargs): super(Material, self).__init__(name, parser, *args, **kwargs) - self.is_non_periodic = kwargs.get("is_non_periodic", False) + self.is_non_periodic = kwargs.get("is_non_periodic", getattr(parser, "is_non_periodic", False)) cell_type = kwargs.get("cell_type", "original") structure_string = kwargs.get("structure_string") @@ -36,6 +37,8 @@ def __init__(self, name, parser, *args, **kwargs): else: basis = self.parser.initial_basis() lattice = self.parser.initial_lattice_vectors() + if self.is_non_periodic and lattice is None: + lattice, basis = box_molecule(basis, [basis, self.parser.final_basis()]) structure_string = lattice_basis_to_poscar(lattice, basis) if kwargs.get("is_final_structure"): @@ -45,6 +48,8 @@ def __init__(self, name, parser, *args, **kwargs): else: basis = self.parser.final_basis() lattice = self.parser.final_lattice_vectors() + if self.is_non_periodic and lattice is None: + lattice, basis = box_molecule(basis, [self.parser.initial_basis(), basis]) structure_string = lattice_basis_to_poscar(lattice, basis) if self.is_non_periodic: @@ -114,11 +119,12 @@ def _serialize(self): "unitCellFormula": self.unitCellFormula, "lattice": self.lattice, "basis": self.basis, + "isNonPeriodic": self.is_non_periodic, "derivedProperties": self.derived_properties, "creator": {"_id": "", "cls": "User", "slug": ""}, "owner": {"_id": "", "cls": "Account", "slug": ""}, "schemaVersion": "0.2.0", - "metadata": {} + "metadata": {}, } def _elemental_ratios(self): diff --git a/express/properties/utils.py b/express/properties/utils.py index d91d4892..4ecdda81 100644 --- a/express/properties/utils.py +++ b/express/properties/utils.py @@ -1,3 +1,7 @@ +from mat3ra.made.tools.convert.utils import calculate_padded_cell_simple_cubic +from mat3ra.made.utils import get_center_of_coordinates + + def eigenvalues(eigenvalues_at_kpoints, kpoint_index=0, spin_index=0): """ Returns eigenvalues for a given kpoint and spin. @@ -27,3 +31,23 @@ def to_array_with_ids(array): list """ return [{"id": index, "value": value} for index, value in enumerate(array)] + + +def box_molecule(selected_basis, parsed_bases): + """ + Returns lattice vectors and the selected basis centered in them: made's simple-cubic padding, + sized to hold every basis in `parsed_bases` so an optimization moves atoms inside a fixed box + rather than resizing it. Printed coordinates straddle the origin, and atoms left outside the box + read as extra fragments and corrupt the InChI. + """ + coordinates = [coordinate["value"] for coordinate in selected_basis["coordinates"]] + cells = [ + calculate_padded_cell_simple_cubic([point["value"] for point in other["coordinates"]]) for other in parsed_bases + ] + edge = max(vectors[0][0] for vectors in cells) + center = get_center_of_coordinates(coordinates) + centered = [[x - center[axis] + edge / 2 for axis, x in enumerate(coordinate)] for coordinate in coordinates] + return ( + {"vectors": {"a": [edge, 0.0, 0.0], "b": [0.0, edge, 0.0], "c": [0.0, 0.0, edge], "alat": 1}}, + dict(selected_basis, coordinates=to_array_with_ids(centered)), + ) diff --git a/tests/fixtures/nwchem/references.py b/tests/fixtures/nwchem/references.py index abb10958..20f38c8c 100644 --- a/tests/fixtures/nwchem/references.py +++ b/tests/fixtures/nwchem/references.py @@ -39,12 +39,22 @@ "nuclear_repulsion": {"name": "nuclear_repulsion", "value": 250.20815670232923}, } +# test-001 declares `units au`, so these are its printed coordinates converted to angstrom. BASIS = { "units": "angstrom", "elements": [{"id": 0, "value": "O"}, {"id": 1, "value": "H"}, {"id": 2, "value": "H"}], "coordinates": [ - {"id": 0, "value": [0.00000000, 0.00000000, 0.22143053]}, - {"id": 1, "value": [0.00000000, 1.43042809, -0.88572213]}, - {"id": 2, "value": [0.00000000, -1.43042809, -0.88572213]}, + {"id": 0, "value": [0.0, 0.0, 0.11717600]}, + {"id": 1, "value": [0.0, 0.75695001, -0.46870401]}, + {"id": 2, "value": [0.0, -0.75695001, -0.46870401]}, ], } + +# What express serializes for test-002: made's cubic padding sized to the larger of its two blocks, +# with the final basis centered in it. The edge therefore pins the initial block too. +FINAL_CELL_EDGE_MULTISTEP = 5.721712 +FINAL_CRYSTAL_COORDINATES_MULTISTEP = [ + [0.5, 0.5, 0.569589977], + [0.366715633, 0.5, 0.465205011], + [0.633284367, 0.5, 0.465205011], +] diff --git a/tests/integration/parsers/apps/nwchem/test_parser.py b/tests/integration/parsers/apps/nwchem/test_parser.py index a0fbc66b..85c63c4f 100644 --- a/tests/integration/parsers/apps/nwchem/test_parser.py +++ b/tests/integration/parsers/apps/nwchem/test_parser.py @@ -47,6 +47,10 @@ def test_nwchem_lumo_energy_multistep(self): self.assertAlmostEqual(lumo_energy, LUMO_ENERGY_MULTISTEP, places=2) self.assertNotAlmostEqual(lumo_energy, LUMO_ENERGY_MULTISTEP_INITIAL_GUESS, places=2) + def test_nwchem_structures_of_single_point(self): + self.assertDeepAlmostEqual(self.parser.initial_basis(), BASIS, places=6) + self.assertDeepAlmostEqual(self.parser.final_basis(), BASIS, places=6) + def test_nwchem_total_energy_contributions(self): self.assertDeepAlmostEqual(self.parser.total_energy_contributions(), TOTAL_ENERGY_CONTRIBUTION, places=2) diff --git a/tests/integration/properties/test_material.py b/tests/integration/properties/test_material.py index a3427fcb..3512b20c 100644 --- a/tests/integration/properties/test_material.py +++ b/tests/integration/properties/test_material.py @@ -3,9 +3,11 @@ from typing import Dict, List from express.parsers.apps.espresso.parser import EspressoParser +from express.parsers.apps.nwchem.parser import NwchemParser from express.parsers.apps.vasp.parser import VaspParser from express.properties.material import Material from tests.fixtures.data import SI as data +from tests.fixtures.nwchem.references import FINAL_CELL_EDGE_MULTISTEP, FINAL_CRYSTAL_COORDINATES_MULTISTEP from tests.integration import IntegrationTestBase @@ -24,6 +26,10 @@ def vasp_parser(self): def espresso_parser(self): return EspressoParser(work_dir=self.workDir, stdout_file=self.stdoutFile) + @property + def nwchem_parser(self): + return NwchemParser(work_dir=self.workDir, stdout_file=self.stdoutFile) + @property def structure_string(self): manifest = self.getManifest() @@ -46,6 +52,9 @@ def assertJsonEqual(self, material: Material) -> bool: derived_props = self.filter_derived_props(material.is_non_periodic) json = deepcopy(data) json["derivedProperties"] = derived_props + # Same reason derivedProperties is adapted: the shared SI fixture is periodic, and each + # test picks the picture it wants. + json["isNonPeriodic"] = material.is_non_periodic self.assertDeepAlmostEqual(material.serialize_and_validate(), json, places=2) return True @@ -93,3 +102,15 @@ def test_material_from_structure(self): def test_material_serialize_and_validate(self): material = Material("material", self.vasp_parser, is_initial_structure=True, is_non_periodic=True) self.assertJsonEqual(material) + + def test_nwchem_material_of_a_relaxed_molecule(self): + # Constructed WITHOUT is_non_periodic on purpose: rupy never passes it. + initial = Material("material", self.nwchem_parser, is_initial_structure=True).serialize_and_validate() + material = Material("material", self.nwchem_parser, is_final_structure=True).serialize_and_validate() + coordinates = [c["value"] for c in material["basis"]["coordinates"]] + self.assertEqual(initial["lattice"], material["lattice"]) + self.assertEqual([material["isNonPeriodic"], material["lattice"]["type"]], [True, "CUB"]) + self.assertAlmostEqual(material["lattice"]["a"], FINAL_CELL_EDGE_MULTISTEP, places=6) + self.assertDeepAlmostEqual(coordinates, FINAL_CRYSTAL_COORDINATES_MULTISTEP, places=6) + derived = {p["name"] for p in material["derivedProperties"]} + self.assertEqual(derived & {"inchi", "inchi_key", "volume", "density"}, {"inchi", "inchi_key"}) diff --git a/tests/manifest.yaml b/tests/manifest.yaml index aa9c6001..fad9041c 100644 --- a/tests/manifest.yaml +++ b/tests/manifest.yaml @@ -30,6 +30,10 @@ test_nwchem_total_energy_contributions: workDir: fixtures/nwchem/test-001 stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log +test_nwchem_structures_of_single_point: + workDir: fixtures/nwchem/test-001 + stdoutFile: fixtures/nwchem/test-001/nwchem-total-energy.log + test_nwchem_zero_point_energy: workDir: fixtures/nwchem/test-002 stdoutFile: fixtures/nwchem/test-002/nwchem-frequency.log @@ -256,3 +260,7 @@ test_espresso_hubbard_v: test_espresso_hubbard_v_nn: workDir: fixtures/espresso/v7_2/test-009 stdoutFile: fixtures/espresso/v7_2/test-009/HUBBARD.dat + +test_nwchem_material_of_a_relaxed_molecule: + workDir: fixtures/nwchem/test-002 + stdoutFile: fixtures/nwchem/test-002/nwchem-frequency.log diff --git a/tests/unit/parsers/test_nwchem_txt_parser.py b/tests/unit/parsers/test_nwchem_txt_parser.py index a3256d54..9c3d8c86 100644 --- a/tests/unit/parsers/test_nwchem_txt_parser.py +++ b/tests/unit/parsers/test_nwchem_txt_parser.py @@ -1,3 +1,4 @@ +import math import unittest from express.parsers.apps.nwchem.formats.txt import NwchemTXTParser @@ -58,3 +59,69 @@ def test_reads_last_step_of_each_channel(self): def test_returns_nothing_without_an_orbital_analysis_section(self): self.assertEqual(self.parser.eigenvalues_at_vectors(" Total DFT energy = -76.4\n"), []) + +# NOTE: use fixtures in the future instead of inline strings below; tentatively, mat3ra-fixtures +def geometry_block(units, scale, *rows): + return "\n".join( + [ + f" Output coordinates in {units} (scale by {scale} to convert to a.u.)", + "", + " No. Tag Charge X Y Z", + " ---- ---------------- ---------- -------------- -------------- --------------", + *rows, + "", + " Atomic Mass ", + " O 15.994910", + "", + ] + ) + + +ANGSTROM = ("angstroms", "1.889725989") +ATOMIC_UNITS = ("a.u.", "1.000000000") +ANGSTROM_ROWS = [ + " 1 O 8.0000 0.00000000 0.00000000 1.00000000", + " 2 H 1.0000 0.00000000 0.00000000 -1.00000000", +] + + +class NwchemTXTParserGeometryTest(unittest.TestCase): + def setUp(self): + self.parser = NwchemTXTParser(work_dir=".") + + def coordinates(self, text, index=-1): + return [[round(v, 6) for v in c["value"]] for c in self.parser.basis(text, index)["coordinates"]] + + def test_selects_the_block_by_index(self): + text = geometry_block(*ANGSTROM, *ANGSTROM_ROWS) + geometry_block( + *ANGSTROM, + " 1 O 8.0000 0.00000000 0.00000000 0.50000000", + " 2 H 1.0000 0.00000000 0.00000000 -0.50000000", + ) + self.assertEqual([self.coordinates(text, 0)[0][2], self.coordinates(text)[0][2]], [1.0, 0.5]) + + def test_bond_length_of_an_atomic_units_block_is_physical(self): + # test-001's rows verbatim. Read as angstrom they give an O-H of 1.81 A, not 0.9572 A. + text = geometry_block( + *ATOMIC_UNITS, + " 1 O 8.0000 0.00000000 0.00000000 0.22143053", + " 2 H 1.0000 0.00000000 1.43042809 -0.88572213", + " 3 H 1.0000 0.00000000 -1.43042809 -0.88572213", + ) + coordinates = self.coordinates(text) + self.assertAlmostEqual(math.dist(coordinates[0], coordinates[1]), 0.9572, places=4) + self.assertAlmostEqual(math.dist(coordinates[0], coordinates[2]), 0.9572, places=4) + + def test_refuses_blocks_that_disagree_on_their_atoms(self): + # A log truncated mid-table parses as a shorter molecule, which rupy would publish. + text = geometry_block(*ANGSTROM, *ANGSTROM_ROWS) + geometry_block(*ANGSTROM, ANGSTROM_ROWS[0]) + self.assertIsNone(self.parser.basis(text, -1)) + + def test_returns_nothing_without_a_geometry_block(self): + self.assertIsNone(self.parser.basis(" Total DFT energy = -76.4\n", -1)) + + def test_ignores_a_numeric_row_that_happens_to_fit_the_shape(self): + # A row of bare numbers satisfies the column shape but has no element symbol. It must be + # skipped, not raise: rupy swallows the exception and final_structure vanishes silently. + text = geometry_block(*ANGSTROM, *ANGSTROM_ROWS, " 3 1.0000 2.0000 3.0000 4.0000 5.0000") + self.assertEqual([e["value"] for e in self.parser.basis(text, -1)["elements"]], ["O", "H"]) diff --git a/tests/unit/properties/test_utils.py b/tests/unit/properties/test_utils.py new file mode 100644 index 00000000..c757610a --- /dev/null +++ b/tests/unit/properties/test_utils.py @@ -0,0 +1,45 @@ +import unittest + +from express.parsers.mixins.ionic import IonicDataMixin +from express.properties.material import Material +from express.properties.utils import box_molecule, to_array_with_ids + + +def basis(*heights): + return { + "units": "angstrom", + "elements": to_array_with_ids(["H"] * len(heights)), + "coordinates": to_array_with_ids([[0.0, 0.0, height] for height in heights]), + } + + +COMPACT, EXPANDED = basis(0.2, -0.2), basis(2.5, -2.5) +EXPANDED_EDGE = 15.0 +STRUCTURE_KWARGS = ("is_initial_structure", "is_final_structure") + + +class ExpandingRelaxationParser(IonicDataMixin): + is_non_periodic = True + + def initial_basis(self): + return COMPACT + + def final_basis(self): + return EXPANDED + + +class BoxMoleculeTest(unittest.TestCase): + def test_shared_cell_holds_a_relaxation_that_expands(self): + # No committed fixture relaxes outward, and that is the direction that leaves atoms outside + # the box, where they read as extra fragments and corrupt the InChI. + for selected_basis in (COMPACT, EXPANDED): + lattice, centered = box_molecule(selected_basis, [COMPACT, EXPANDED]) + self.assertEqual(lattice["vectors"]["a"][0], EXPANDED_EDGE) + for coordinate in centered["coordinates"]: + self.assertTrue(all(0.0 <= x <= EXPANDED_EDGE for x in coordinate["value"])) + + def test_material_hands_over_every_parsed_basis(self): + parser = ExpandingRelaxationParser() + lattices = [Material("material", parser, **{kwarg: True}).lattice for kwarg in STRUCTURE_KWARGS] + self.assertEqual(lattices[0], lattices[1]) + self.assertAlmostEqual(lattices[0]["a"], EXPANDED_EDGE, places=6)