Improved the xTB parser (#856)

Author: LeenFahoum

arc/constants.pxd

@@ -1 +1 @@
-cdef double pi, Na, kB, R, h, hbar, c, e, m_e, m_p, m_n, amu, a0, E_h, F, E_h_kJmol, bohr_to_angstrom
+cdef double pi, Na, kB, R, h, hbar, c, e, m_e, m_p, m_n, amu, a0, bohr_to_angstrom, E_h, F, E_h_kJmol

arc/constants.py

@@ -41,6 +41,7 @@
 
 #: The Bohr radius :math:`a_0` in :math:`\mathrm{m}`
 a0 = 5.2917721092e-11
+bohr_to_angstrom = a0 * 1e10
 
 #: The atomic mass unit in :math:`\mathrm{kg}`
 amu = 1.660538921e-27
@@ -78,8 +79,6 @@
 #: Vacuum permittivity
 epsilon_0 = 8.8541878128
 
-bohr_to_angstrom = 0.529177
-
 # Cython does not automatically place module-level variables into the module
 # symbol table when in compiled mode, so we must do this manually so that we
 # can use the constants from both Cython and regular Python code

arc/parser/adapters/xtb.py

@@ -11,7 +11,7 @@
 import pandas as pd
 
 from arc.common import is_str_float
-from arc.constants import E_h_kJmol
+from arc.constants import E_h_kJmol, bohr_to_angstrom
 from arc.species.converter import str_to_xyz, xyz_from_data, logger
 from arc.parser.adapter import ESSAdapter
 from arc.parser.factory import register_ess_adapter
@@ -35,13 +35,22 @@ def logfile_contains_errors(self) -> Optional[str]:
         Returns: Optional[str]
             ``None`` if the log file is free of errors, otherwise the error is returned as a string.
         """
-        # Not implemented for xTB.
+        with open(self.log_file_path, 'r') as f:
+            for line in f:
+                if 'abnormal termination' in line.lower():
+                    return line.strip()
+                if '[ERROR]' in line or '#ERROR' in line:
+                    return line.strip()
         return None
 
     def parse_geometry(self) -> Optional[Dict[str, tuple]]:
         """
         Parse the xyz geometry from an ESS log file.
 
+        Supports both Turbomol ``$coord`` (Bohr) and Molfile V2000 (Angstrom) formats.
+        If the file contains multiple geometry blocks (e.g., from multiple optimization
+        cycles), only the last one is returned.
+
         Returns: Optional[Dict[str, tuple]]
             The cartesian geometry.
         """
@@ -58,6 +67,8 @@ def parse_geometry(self) -> Optional[Dict[str, tuple]]:
                 final_structure = True
                 continue
             if final_structure and '$coord' in line:
+                # Reset on each new $coord block so we keep only the last one
+                coords, symbols = list(), list()
                 in_coord_block = True
                 continue
             if final_structure and 'V2000' in line and not in_coord_block:
@@ -66,16 +77,18 @@ def parse_geometry(self) -> Optional[Dict[str, tuple]]:
                 if parts and parts[0].isdigit():
                     atom_count = int(parts[0])
                     molfile_line_counter = 0
+                # Reset on each new V2000 block
+                coords, symbols = list(), list()
                 continue
 
-            # Parse $coord format
+            # Parse $coord format (Turbomole $coord coordinates are in Bohr, convert to Angstrom)
             if in_coord_block:
                 if '$' in line or 'end' in line.lower() or len(line.split()) < 4:
                     in_coord_block = False
                     continue
                 parts = line.split()
                 try:
-                    x, y, z = map(float, parts[:3])
+                    x, y, z = (float(v) * bohr_to_angstrom for v in parts[:3])
                     symbol = parts[3].capitalize() if len(parts[3]) == 1 else parts[3][0].upper() + parts[3][1:].lower()
                     coords.append([x, y, z])
                     symbols.append(symbol)
@@ -105,30 +118,34 @@ def parse_frequencies(self) -> Optional[np.ndarray]:
         """
         Parse the frequencies from a freq job output file.
 
+        xTB prints frequencies twice (once after the Hessian and once in the
+        Frequency Printout section). This method reads ALL eigval blocks and
+        returns only the last complete one to ensure we get the final values.
+
         Returns: Optional[np.ndarray]
             The parsed frequencies (in cm^-1).
         """
-        freqs = list()
+        # Collect all eigval blocks; use the last one
+        all_blocks = list()
+        current_block = list()
         lines = _get_lines_from_file(self.log_file_path)
-        read_output = False
 
         for line in lines:
-            if read_output:
-                if 'eigval :' in line:
-                    splits = line.split()
-                    for split in splits[2:]:
-                        try:
-                            freq = float(split)
-                            if freq != 0.0:
-                                freqs.append(freq)
-                        except ValueError:
-                            continue
-                elif line.strip() == "" or "projected vibrational frequencies" in line.lower():
-                    continue
-                else:
-                    break
-            if 'vibrational frequencies' in line.lower():
-                read_output = True
+            if 'eigval :' in line:
+                splits = line.split()
+                for split in splits[2:]:
+                    try:
+                        current_block.append(float(split))
+                    except ValueError:
+                        continue
+            elif current_block:
+                # End of an eigval run
+                all_blocks.append(current_block)
+                current_block = list()
+        if current_block:
+            all_blocks.append(current_block)
+
+        freqs = [f for f in all_blocks[-1] if f != 0.0] if all_blocks else list()
 
         # Fallback: try vibspectrum file if no frequencies found in output
         if not freqs:
@@ -232,28 +249,32 @@ def parse_e_elect(self) -> Optional[float]:
         """
         Parse the electronic energy from an sp job output file.
 
+        Looks for ``:: total energy ... Eh`` (SUMMARY block) or
+        ``| TOTAL ENERGY ... Eh |`` (TOTAL section).
+        Avoids false matches against ``total energy gain`` (optimization deltas).
+
         Returns: Optional[float]
             The electronic energy in kJ/mol.
         """
+        import re
         lines = _get_lines_from_file(self.log_file_path)
         energy = None
-        for line in reversed(lines):
-            if 'total energy' in line.lower():
-                try:
-                    energy = float(line.split()[3].strip())
-                    break
-                except (ValueError, IndexError):
+        # Iterate forward and keep the LAST hit (final result)
+        for line in lines:
+            stripped = line.strip()
+            if stripped.startswith(':: total energy') or 'TOTAL ENERGY' in line:
+                m = re.search(r'(-?\d+\.\d+)\s+Eh', line)
+                if m:
+                    energy = float(m.group(1))
+        if energy is None:
+            # Fallback: 'final energy' lines (rare)
+            for line in reversed(lines):
+                if 'final energy' in line.lower():
                     try:
-                        energy = float(line.split()[-1].strip())
+                        energy = float(line.split()[-1])
                         break
                     except (ValueError, IndexError):
                         continue
-            if 'final energy' in line.lower():
-                try:
-                    energy = float(line.split()[-1])
-                    break
-                except (ValueError, IndexError):
-                    continue
         if energy is not None:
             return energy * E_h_kJmol
         return None
@@ -265,16 +286,16 @@ def parse_zpe_correction(self) -> Optional[float]:
         Returns: Optional[float]
             The calculated zero point energy in kJ/mol.
         """
+        import re
         zpe = None
         with open(self.log_file_path, 'r') as f:
             for line in f:
                 if 'zero-point vibrational energy' in line.lower() or 'zero point energy' in line.lower():
                     #          :: zero point energy           0.056690417480 Eh   ::
-                    try:
-                        zpe = float(line.split()[-3])
+                    m = re.search(r'(\d+\.\d+(?:[eE][+-]?\d+)?)\s+Eh', line)
+                    if m:
+                        zpe = float(m.group(1))
                         break
-                    except (ValueError, IndexError):
-                        continue
         if zpe is not None:
             return zpe * E_h_kJmol
         return None
@@ -323,8 +344,12 @@ def parse_1d_scan_energies(self) -> Tuple[Optional[List[float]], Optional[List[f
             logger.warning(f'No valid scan points found in xTB scan log file {scan_path}.')
             return None, None
 
-        # Angles: evenly spaced 0 to 360 inclusive
-        angles = [i * 360.0 / (n_points + 1) for i in range(n_points + 1)]
+        # Angles: evenly spaced from 0 deg with one angle per energy point.
+        # For 44 energies, the dihedral was scanned in steps of 360/45 = 8 deg.
+        # We return n_points angles (matching the energies length): [0, 8, 16, ..., 344].
+        # Note: callers expecting (n_points+1) angles for n_points energies should
+        # add the closing 360 deg themselves.
+        angles = [i * 360.0 / (n_points + 1) for i in range(n_points)]
 
         return rel_energies, angles
 

arc/parser/parser_test.py

@@ -480,18 +480,18 @@ def test_parse_geometry_2(self):
         self.assertIsInstance(xyz_1, dict)
         self.assertEqual(len(xyz_1['symbols']), 53)
 
-        path_2 = os.path.join(ARC_TESTING_PATH, 'opt', 'xtb_opt_1.out')  # Turbomol format
+        path_2 = os.path.join(ARC_TESTING_PATH, 'opt', 'xtb_opt_1.out')  # Turbomol format ($coord is in Bohr)
         xyz_2 = parser.parse_geometry(log_file_path=path_2)
         expected_xyz_2 = {'symbols': ('C', 'C', 'C', 'O', 'H', 'H', 'H', 'H'),
                           'isotopes': (12, 12, 12, 16, 1, 1, 1, 1),
-                          'coords': ((-2.1908361683609, -0.0875055545944093, -0.712358508847116),
-                                     (-0.170663821576948, -0.721009080780027, 0.6339514359292),
-                                     (2.27150340995404, 0.539200343733434, 0.305747355748416),
-                                     (4.16185183882216, 0.0456870795397582, 1.46958543963615),
-                                     (-2.10906550336131, 1.39281538085596, -2.11159862590894),
-                                     (-3.98822887666127, -1.01150815322938, -0.474092464462499),
-                                     (-0.230883305078883, -2.20128348308961, 2.04098984698359),
-                                     (2.25632242626311, 2.04360346756428, -1.15222446018155))}
+                          'coords': ((-1.1593401110647, -0.0463059268636, -0.3769637386362),
+                                     (-0.0903113691106, -0.3815414223399, 0.3354725190107),
+                                     (1.2020273599692, 0.2853324202958, 0.1617944684729),
+                                     (2.2023562705124, 0.0241765516896, 0.7776708141903),
+                                     (-1.1160689558722, 0.7370458647952, -1.1174094260626),
+                                     (-2.110478992265, -0.5352668500015, -0.2508788280669),
+                                     (-0.1221781347317, -1.1648685897309, 1.0800448842572),
+                                     (1.1939939325626, 1.0814279521553, -0.6097306831655))}
         self.assertTrue(almost_equal_coords(xyz_2, expected_xyz_2))
 
         path_3 = os.path.join(ARC_TESTING_PATH, 'opt', 'xtb_opt_2.out')  # SDF format
@@ -732,10 +732,12 @@ def test_parse_1d_scan_energies(self):
                                       8.397974544321187, 7.071194090611243, 5.214457982623571, 3.4362612986915337,
                                       1.958199294316728, 0.8766536693692615, 0.22504856466548517,
                                       0.0004629911018128041])
+        # Length-matched: 44 energies -> 44 angles (step = 360/(n+1) = 8 deg)
+        self.assertEqual(len(angles_3), len(energies_3))
         self.assertEqual(angles_3, [0.0, 8.0, 16.0, 24.0, 32.0, 40.0, 48.0, 56.0, 64.0, 72.0, 80.0, 88.0, 96.0,
                                     104.0, 112.0, 120.0, 128.0, 136.0, 144.0, 152.0, 160.0, 168.0, 176.0, 184.0, 192.0,
                                     200.0, 208.0, 216.0, 224.0, 232.0, 240.0, 248.0, 256.0, 264.0, 272.0, 280.0, 288.0,
-                                    296.0, 304.0, 312.0, 320.0, 328.0, 336.0, 344.0, 352.0])
+                                    296.0, 304.0, 312.0, 320.0, 328.0, 336.0, 344.0])
 
         path_4 = os.path.join(ARC_TESTING_PATH, 'rotor_scans', 'orca', 'cc.txt')
         energies_4, angles_4 = parser.parse_1d_scan_energies(log_file_path=path_4)
@@ -1000,6 +1002,81 @@ def test_parse_active_space(self):
                                                               xyz="""N      0.0   0.0    0.0"""))
         self.assertEqual(active, {'e_o': (5, 4), 'occ': [3, 1, 1, 0, 1, 0, 0, 0], 'closed': [1, 0, 0, 0, 0, 0, 0, 0]})
 
+    # ----- xTB-specific improvements -----
+
+    def test_xtb_logfile_contains_errors(self):
+        """xTB error detection should return None for clean files."""
+        from arc.parser.adapters.xtb import XTBParser
+        clean_path = os.path.join(ARC_TESTING_PATH, 'sp', 'NCC_xTB.out')
+        adapter = XTBParser(log_file_path=clean_path)
+        self.assertIsNone(adapter.logfile_contains_errors())
+
+    def test_xtb_parse_e_elect_avoids_false_match(self):
+        """parse_e_elect should not match 'total energy gain' lines from optimization deltas.
+
+        The CO2_xtb.out file contains both 'total energy gain : -0.1724555 Eh' (delta)
+        and ':: total energy -10.308452243026 Eh' (final). The parser must return the
+        final value, not the delta.
+        """
+        from arc.constants import E_h_kJmol
+        path = os.path.join(ARC_TESTING_PATH, 'freq', 'CO2_xtb.out')
+        e_elect = parser.parse_e_elect(path)
+        self.assertAlmostEqual(e_elect, -10.308452243026 * E_h_kJmol, places=3)
+        # Make sure we did NOT pick up the delta value
+        delta_value = -0.1724555 * E_h_kJmol  # ~-452 kJ/mol
+        self.assertNotAlmostEqual(e_elect, delta_value, places=1)
+
+    def test_xtb_parse_frequencies_takes_last_block(self):
+        """parse_frequencies should return the LAST eigval block.
+
+        TS_NH2+N2H3_xtb.out prints frequencies twice. Both blocks should be identical;
+        if the parser stops after the first block, it would still get the right values,
+        but using the last block is more robust against truncated/duplicate blocks.
+        """
+        path = os.path.join(ARC_TESTING_PATH, 'freq', 'TS_NH2+N2H3_xtb.out')
+        freqs = parser.parse_frequencies(log_file_path=path)
+        self.assertIsNotNone(freqs)
+        # Imaginary frequency should be the first non-zero one
+        self.assertAlmostEqual(freqs[0], -781.89, places=1)
+        # Last frequency
+        self.assertAlmostEqual(freqs[-1], 3467.59, places=1)
+        # 24 modes total (3*8 - 6 = 18 + 6 zero = 24, but only 18 non-zero in TS)
+        self.assertEqual(len(freqs), 18)
+
+    def test_xtb_parse_zpe_correction_regex(self):
+        """parse_zpe_correction uses regex to find ZPE in scientific notation safely."""
+        from arc.constants import E_h_kJmol
+        path = os.path.join(ARC_TESTING_PATH, 'freq', 'TS_NH2+N2H3_xtb.out')
+        zpe = parser.parse_zpe_correction(log_file_path=path)
+        self.assertIsNotNone(zpe)
+        self.assertAlmostEqual(zpe, 0.056690417480 * E_h_kJmol, places=3)
+
+    def test_xtb_parse_1d_scan_energies_length_match(self):
+        """Energies and angles returned by parse_1d_scan_energies must have matching lengths."""
+        path = os.path.join(ARC_TESTING_PATH, 'rotor_scans', 'xtb_1', 'output.out')
+        energies, angles = parser.parse_1d_scan_energies(log_file_path=path)
+        self.assertIsNotNone(energies)
+        self.assertIsNotNone(angles)
+        self.assertEqual(len(energies), len(angles))
+        # First angle = 0, all angles strictly increasing
+        self.assertEqual(angles[0], 0.0)
+        for i in range(1, len(angles)):
+            self.assertGreater(angles[i], angles[i - 1])
+        # All angles should be in [0, 360)
+        self.assertLess(angles[-1], 360.0)
+
+    def test_xtb_parse_geometry_resets_on_new_block(self):
+        """parse_geometry should not accumulate atoms across multiple geometry blocks.
+
+        For CO2_xtb.out (Turbomol $coord with 3 atoms), the parser should return
+        exactly 3 atoms, not 6 or more even if the file mentions multiple blocks.
+        """
+        path = os.path.join(ARC_TESTING_PATH, 'freq', 'CO2_xtb.out')
+        xyz = parser.parse_geometry(log_file_path=path)
+        self.assertIsNotNone(xyz)
+        self.assertEqual(len(xyz['symbols']), 3)
+        self.assertEqual(xyz['symbols'], ('O', 'C', 'O'))
+
 
     def test_parse_opt_steps(self):
         """Test parsing the number of optimization steps from various ESS output files."""