Add multipole parsing for Q-Chem IO (#3490)

* First attempt at multipole moment parsing * Multipole parsing appears to work * Minor docstring change * Linting * Requested changes * more refactor --------- Co-authored-by: Janosh Riebesell <[email protected]>
materialsproject · Dec 5, 2023 · 8a4e86d · 8a4e86d
1 parent a48d541
commit 8a4e86d
Show file tree

Hide file tree

Showing 4 changed files with 101 additions and 3 deletions.
diff --git a/pymatgen/io/qchem/outputs.py b/pymatgen/io/qchem/outputs.py
@@ -948,9 +948,10 @@ def _read_SCF(self):
     def _read_charges_and_dipoles(self):
         """
         Parses Mulliken/ESP/RESP charges.
-        Parses associated dipoles.
+        Parses associated dipole/multipole moments.
         Also parses spins given an unrestricted SCF.
         """
+
         self.data["dipoles"] = {}
         temp_dipole_total = read_pattern(
             self.text, {"key": r"X\s*[\d\-\.]+\s*Y\s*[\d\-\.]+\s*Z\s*[\d\-\.]+\s*Tot\s*([\d\-\.]+)"}
@@ -976,6 +977,81 @@ def _read_charges_and_dipoles(self):
                         dipole[ii][jj] = temp_dipole[ii][jj]
                 self.data["dipoles"]["dipole"] = dipole
 
+        self.data["multipoles"] = dict()
+
+        quad_mom_pat = (
+            r"\s*Quadrupole Moments \(Debye\-Ang\)\s+XX\s+([\-\.0-9]+)\s+XY\s+([\-\.0-9]+)\s+YY"
+            r"\s+([\-\.0-9]+)\s+XZ\s+([\-\.0-9]+)\s+YZ\s+([\-\.0-9]+)\s+ZZ\s+([\-\.0-9]+)"
+        )
+        temp_quadrupole_moment = read_pattern(self.text, {"key": quad_mom_pat}).get("key")
+        if temp_quadrupole_moment is not None:
+            keys = ("XX", "XY", "YY", "XZ", "YZ", "ZZ")
+            if len(temp_quadrupole_moment) == 1:
+                self.data["multipoles"]["quadrupole"] = {
+                    key: float(temp_quadrupole_moment[0][idx]) for idx, key in enumerate(keys)
+                }
+            else:
+                self.data["multipoles"]["quadrupole"] = list()
+                for qpole in temp_quadrupole_moment:
+                    self.data["multipoles"]["quadrupole"].append(
+                        {key: float(qpole[idx]) for idx, key in enumerate(keys)}
+                    )
+
+        octo_mom_pat = (
+            r"\s*Octopole Moments \(Debye\-Ang\^2\)\s+XXX\s+([\-\.0-9]+)\s+XXY\s+([\-\.0-9]+)"
+            r"\s+XYY\s+([\-\.0-9]+)\s+YYY\s+([\-\.0-9]+)\s+XXZ\s+([\-\.0-9]+)\s+XYZ\s+([\-\.0-9]+)"
+            r"\s+YYZ\s+([\-\.0-9]+)\s+XZZ\s+([\-\.0-9]+)\s+YZZ\s+([\-\.0-9]+)\s+ZZZ\s+([\-\.0-9]+)"
+        )
+        temp_octopole_moment = read_pattern(self.text, {"key": octo_mom_pat}).get("key")
+        if temp_octopole_moment is not None:
+            keys = ("XXX", "XXY", "XYY", "YYY", "XXZ", "XYZ", "YYZ", "XZZ", "YZZ", "ZZZ")
+            if len(temp_octopole_moment) == 1:
+                self.data["multipoles"]["octopole"] = {
+                    key: float(temp_octopole_moment[0][idx]) for idx, key in enumerate(keys)
+                }
+            else:
+                self.data["multipoles"]["octopole"] = list()
+                for opole in temp_octopole_moment:
+                    self.data["multipoles"]["octopole"].append({key: float(opole[idx]) for idx, key in enumerate(keys)})
+
+        hexadeca_mom_pat = (
+            r"\s*Hexadecapole Moments \(Debye\-Ang\^3\)\s+XXXX\s+([\-\.0-9]+)\s+XXXY\s+([\-\.0-9]+)"
+            r"\s+XXYY\s+([\-\.0-9]+)\s+XYYY\s+([\-\.0-9]+)\s+YYYY\s+([\-\.0-9]+)\s+XXXZ\s+([\-\.0-9]+)"
+            r"\s+XXYZ\s+([\-\.0-9]+)\s+XYYZ\s+([\-\.0-9]+)\s+YYYZ\s+([\-\.0-9]+)\s+XXZZ\s+([\-\.0-9]+)"
+            r"\s+XYZZ\s+([\-\.0-9]+)\s+YYZZ\s+([\-\.0-9]+)\s+XZZZ\s+([\-\.0-9]+)\s+YZZZ\s+([\-\.0-9]+)"
+            r"\s+ZZZZ\s+([\-\.0-9]+)"
+        )
+        temp_hexadecapole_moment = read_pattern(self.text, {"key": hexadeca_mom_pat}).get("key")
+        if temp_hexadecapole_moment is not None:
+            keys = (
+                "XXXX",
+                "XXXY",
+                "XXYY",
+                "XYYY",
+                "YYYY",
+                "XXXZ",
+                "XXYZ",
+                "XYYZ",
+                "YYYZ",
+                "XXZZ",
+                "XYZZ",
+                "YYZZ",
+                "XZZZ",
+                "YZZZ",
+                "ZZZZ",
+            )
+
+            if len(temp_hexadecapole_moment) == 1:
+                self.data["multipoles"]["hexadecapole"] = {
+                    key: float(temp_hexadecapole_moment[0][idx]) for idx, key in enumerate(keys)
+                }
+            else:
+                self.data["multipoles"]["hexadecapole"] = list()
+                for hpole in temp_hexadecapole_moment:
+                    self.data["multipoles"]["hexadecapole"].append(
+                        {key: float(hpole[idx]) for idx, key in enumerate(keys)}
+                    )
+
         if self.data.get("unrestricted", []):
             header_pattern = (
                 r"\-+\s+Ground-State Mulliken Net Atomic Charges\s+Atom\s+Charge \(a\.u\.\)\s+"

diff --git a/tests/io/qchem/multi_job.json b/tests/io/qchem/multi_job.json
diff --git a/tests/io/qchem/single_job.json b/tests/io/qchem/single_job.json
diff --git a/tests/io/qchem/test_outputs.py b/tests/io/qchem/test_outputs.py
@@ -148,6 +148,7 @@
     "norm_of_stepsize",
     "version",
     "dipoles",
+    "multipoles",
     "gap_info",
 }
 
@@ -310,6 +311,27 @@ def test_all(self):
         for key in property_list:
             self._test_property(key, single_outs, multi_outs)
 
+    def test_multipole_parsing(self):
+        sp = QCOutput(f"{TEST_FILES_DIR}/molecules/new_qchem_files/nbo.qout")
+
+        mpoles = sp.data["multipoles"]
+        assert len(mpoles["quadrupole"]) == 6
+        assert len(mpoles["octopole"]) == 10
+        assert len(mpoles["hexadecapole"]) == 15
+        assert mpoles["quadrupole"]["XX"] == -51.3957
+        assert mpoles["quadrupole"]["YZ"] == 3.5356
+        assert mpoles["octopole"]["XYY"] == -15.0294
+        assert mpoles["octopole"]["XZZ"] == -14.9756
+        assert mpoles["hexadecapole"]["YYYY"] == -326.317
+        assert mpoles["hexadecapole"]["XYZZ"] == 58.0584
+
+        opt = QCOutput(f"{TEST_FILES_DIR}/molecules/new_qchem_files/ts.out")
+        mpoles = opt.data["multipoles"]
+
+        assert len(mpoles["quadrupole"]) == 5
+        assert len(mpoles["octopole"]) == 5
+        assert len(mpoles["hexadecapole"]) == 5
+
     @unittest.skipIf((openbabel is None), "OpenBabel not installed.")
     def test_structural_change(self):
         t1 = Molecule.from_file(f"{TEST_FILES_DIR}/molecules/structural_change/t1.xyz")