SOF-8005: parse NWChem initial and final geometry

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

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

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).

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<units>\S+) \(scale by\s+(?P<scale>[\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<rows>(?:{})+)".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<element>{})\S*".format(ELEMENT_REGEX),
+        double=DOUBLE_REGEX,
+        x=r"(?P<x>{})".format(DOUBLE_REGEX),
+        y=r"(?P<y>{})".format(DOUBLE_REGEX),
+        z=r"(?P<z>{})".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 (?:(?P<spin>Alpha|Beta) )?Molecular Orbital Analysis")
 VECTOR_REGEX = re.compile(

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):

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)),
+    )

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],
+]

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)
 

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"})

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