If constrains are set in the refinement step now the same constrains are set in the postprocess via the usage of xtb with orca as a driver.

src/mindlessgen/generator/main.py

@@ -460,7 +460,7 @@ def setup_engines(
                 raise ImportError("orca not found.")
         except ImportError as e:
             raise ImportError("orca not found.") from e
-        return ORCA(path, cfg.orca)
+        return ORCA(path, cfg.orca, cfg.xtb)
     elif engine_type == "turbomole":
         try:
             jobex_path = jobex_path_func(cfg.turbomole.jobex_path)

src/mindlessgen/qm/orca.py

@@ -2,22 +2,26 @@
 This module handles all ORCA-related functionality.
 """
 
+from collections import defaultdict
 from pathlib import Path
 import shutil
 import subprocess as sp
 from tempfile import TemporaryDirectory
 
 from ..molecules import Molecule
-from ..prog import ORCAConfig
+from ..prog import DistanceConstraint, ORCAConfig, XTBConfig
 from .base import QMMethod
+from .xtb import get_xtb_path
 
 
 class ORCA(QMMethod):
     """
     This class handles all interaction with the ORCA external dependency.
     """
 
-    def __init__(self, path: str | Path, orcacfg: ORCAConfig) -> None:
+    def __init__(
+        self, path: str | Path, orcacfg: ORCAConfig, xtb_config: XTBConfig | None = None
+    ) -> None:
         """
         Initialize the ORCA class.
         """
@@ -28,6 +32,7 @@ def __init__(self, path: str | Path, orcacfg: ORCAConfig) -> None:
         else:
             raise TypeError("orca_path should be a string or a Path object.")
         self.cfg = orcacfg
+        self.xtb_cfg = xtb_config
         # must be explicitly initialized in current parallelization implementation
         # as accessing parent class variables might not be possible
         self.tmp_dir = self.__class__.get_temporary_directory()
@@ -51,11 +56,22 @@ def optimize(
             # NOTE: "prefix" and "dir" are valid keyword arguments for TemporaryDirectory
             temp_path = Path(temp_dir).resolve()
             # write the molecule to a temporary file
-            molecule.write_xyz_to_file(temp_path / "molecule.xyz")
+            xyz_filename = "molecule.xyz"
+            molecule.write_xyz_to_file(temp_path / xyz_filename)
 
             inputname = "orca_opt.inp"
+            use_xtb_driver = self._should_use_xtb_driver()
+            xtb_input = temp_path / "xtb.inp"
+            if use_xtb_driver:
+                self._write_xtb_input(molecule, xtb_input, inputname)
             orca_input = self._gen_input(
-                molecule, "molecule.xyz", ncores, True, max_cycles
+                molecule,
+                xyz_filename,
+                temp_path,
+                ncores,
+                True,
+                max_cycles,
+                use_xtb_driver=use_xtb_driver,
             )
             if verbosity > 1:
                 print("ORCA input file:\n##################")
@@ -65,13 +81,20 @@ def optimize(
                 f.write(orca_input)
 
             # run orca
-            arguments = [
-                inputname,
-            ]
-
-            orca_log_out, orca_log_err, return_code = self._run(
-                temp_path=temp_path, arguments=arguments
-            )
+            if use_xtb_driver:
+                orca_log_out, orca_log_err, return_code = self._run_xtb_driver(
+                    temp_path=temp_path,
+                    geometry_filename=xyz_filename,
+                    xcontrol_name=xtb_input.name,
+                    ncores=ncores,
+                )
+            else:
+                arguments = [
+                    inputname,
+                ]
+                orca_log_out, orca_log_err, return_code = self._run(
+                    temp_path=temp_path, arguments=arguments
+                )
             if verbosity > 2:
                 print(orca_log_out)
             if return_code != 0:
@@ -80,7 +103,14 @@ def optimize(
                 )
 
             # read the optimized molecule from the output file
-            xyzfile = Path(temp_path / inputname).resolve().with_suffix(".xyz")
+            if use_xtb_driver:
+                xyzfile = temp_path / "xtbopt.xyz"
+                if not xyzfile.exists():
+                    raise RuntimeError(
+                        "xTB-driven ORCA optimization did not produce 'xtbopt.xyz'."
+                    )
+            else:
+                xyzfile = Path(temp_path / inputname).resolve().with_suffix(".xyz")
             optimized_molecule = molecule.copy()
             optimized_molecule.read_xyz_from_file(xyzfile)
             return optimized_molecule
@@ -103,10 +133,10 @@ def singlepoint(self, molecule: Molecule, ncores: int, verbosity: int = 1) -> st
 
             # write the input file
             inputname = "orca.inp"
-            orca_input = self._gen_input(molecule, molfile, ncores)
+            orca_input = self._gen_input(molecule, molfile, temp_path, ncores)
             if verbosity > 1:
                 print("ORCA input file:\n##################")
-                print(self._gen_input(molecule, molfile, ncores))
+                print(self._gen_input(molecule, molfile, temp_path, ncores))
                 print("##################")
             with open(temp_path / inputname, "w", encoding="utf8") as f:
                 f.write(orca_input)
@@ -170,13 +200,172 @@ def _run(self, temp_path: Path, arguments: list[str]) -> tuple[str, str, int]:
             orca_log_err = e.stderr.decode("utf8", errors="replace")
             return orca_log_out, orca_log_err, e.returncode
 
+    def _run_xtb_driver(
+        self,
+        temp_path: Path,
+        geometry_filename: str,
+        xcontrol_name: str,
+        ncores: int,
+    ) -> tuple[str, str, int]:
+        """
+        Run the optimization through the xTB external driver when constraints are requested.
+        """
+        xtb_executable = self._get_xtb_executable()
+        arguments = [
+            str(xtb_executable),
+            geometry_filename,
+            "--opt",
+        ]
+        opt_level = getattr(self.cfg, "optlevel", None)
+        if opt_level not in (None, ""):
+            arguments.append(str(opt_level))
+        arguments.extend(["--orca", "-I", xcontrol_name])
+        try:
+            xtb_out = sp.run(
+                arguments,
+                cwd=temp_path,
+                capture_output=True,
+                check=True,
+            )
+            xtb_log_out = xtb_out.stdout.decode("utf8", errors="replace")
+            xtb_log_err = xtb_out.stderr.decode("utf8", errors="replace")
+            return xtb_log_out, xtb_log_err, 0
+        except sp.CalledProcessError as e:
+            xtb_log_out = e.stdout.decode("utf8", errors="replace")
+            xtb_log_err = e.stderr.decode("utf8", errors="replace")
+            return xtb_log_out, xtb_log_err, e.returncode
+
+    def _get_xtb_executable(self) -> Path:
+        """
+        Determine the path to the xTB executable for external ORCA optimizations.
+        """
+        for attr_name in ("xtb_driver_path", "xtb_path"):
+            candidate = getattr(self.cfg, attr_name, None)
+            if candidate:
+                try:
+                    return get_xtb_path(candidate)
+                except ImportError as exc:
+                    raise RuntimeError(
+                        f"xTB executable defined via '{attr_name}' could not be found."
+                    ) from exc
+        try:
+            return get_xtb_path(None)
+        except ImportError as exc:
+            raise RuntimeError(
+                "xTB executable not found. Required for constrained ORCA optimizations."
+            ) from exc
+
+    def _should_use_xtb_driver(self) -> bool:
+        """
+        Determine if the xTB external driver should be used (constraints configured).
+        """
+        return bool(self.xtb_cfg and self.xtb_cfg.distance_constraints)
+
+    def _write_xtb_input(
+        self, molecule: Molecule, xtb_input: Path, input_file: str
+    ) -> None:
+        """
+        Write the xcontrol file containing constraints and ORCA driver info.
+        """
+        if not self.xtb_cfg:
+            raise RuntimeError(
+                "xTB configuration missing but constraints were requested."
+            )
+        constraint_lines = self._prepare_distance_constraint_section(molecule)
+        lines: list[str] = []
+        if constraint_lines:
+            lines.append("$constrain")
+            if self.xtb_cfg.distance_constraint_force_constant is not None:
+                lines.append(
+                    f"  force constant= {self.xtb_cfg.distance_constraint_force_constant}"
+                )
+            lines.extend(constraint_lines)
+            lines.append("$end")
+        lines.append("$external")
+        lines.append(f"  orca input file= {input_file}")
+        lines.append(f"  orca bin= {self.path}")
+        lines.append("$end")
+        xtb_input.write_text("\n".join(lines) + "\n", encoding="utf8")
+
+    def _prepare_distance_constraint_section(self, molecule: Molecule) -> list[str]:
+        """
+        Convert configured distance constraints to xcontrol instructions.
+        """
+        if not self.xtb_cfg or not self.xtb_cfg.distance_constraints:
+            return []
+        element_map: defaultdict[int, list[int]] = defaultdict(list)
+        for idx, atomic_number in enumerate(molecule.ati):
+            element_map[int(atomic_number)].append(idx)
+        constraint_lines: list[str] = []
+        for constraint in self.xtb_cfg.distance_constraints:
+            self._ensure_constraint_atoms_present(element_map, constraint)
+            pairs = self._generate_constraint_pairs(element_map, constraint)
+            if not pairs:
+                raise RuntimeError(
+                    f"No atom pairs found for distance constraint {constraint}."
+                )
+            for first, second in pairs:
+                constraint_lines.append(
+                    f"  distance: {first + 1}, {second + 1}, {constraint.distance:.5f}"
+                )
+        return constraint_lines
+
+    @staticmethod
+    def _generate_constraint_pairs(
+        element_map: dict[int, list[int]], constraint: DistanceConstraint
+    ) -> list[tuple[int, int]]:
+        """
+        Generate index pairs for the provided constraint.
+        """
+        atom_a, atom_b = constraint.atomic_numbers
+        atom_a_idx = atom_a - 1
+        atom_b_idx = atom_b - 1
+        indices_a = element_map.get(atom_a_idx, [])
+        indices_b = element_map.get(atom_b_idx, [])
+
+        if atom_a == atom_b:
+            if len(indices_a) < 2:
+                return []
+            first, second = sorted(indices_a[:2])
+            return [(first, second)]
+
+        if not indices_a or not indices_b:
+            return []
+
+        first, second = indices_a[0], indices_b[0]
+        if first == second:
+            return []
+        if first > second:
+            first, second = second, first
+        return [(first, second)]
+
+    @staticmethod
+    def _ensure_constraint_atoms_present(
+        element_map: dict[int, list[int]], constraint: DistanceConstraint
+    ) -> None:
+        """
+        Validate that the molecule contains enough atoms for the constraint.
+        """
+        for atomic_number, required in constraint.required_counts().items():
+            idx = atomic_number - 1
+            available = len(element_map.get(idx, []))
+            if available < required:
+                symbol = constraint.symbol_for(atomic_number)
+                raise RuntimeError(
+                    f"Distance constraint {constraint} requires at least "
+                    f"{required} atom(s) of {symbol}, but only {available} present."
+                )
+
     def _gen_input(
         self,
         molecule: Molecule,
         xyzfile: str,
+        temp_path: Path,
         ncores: int,
         optimization: bool = False,
         opt_cycles: int | None = None,
+        *,
+        use_xtb_driver: bool = False,
     ) -> str:
         """
         Generate a default input file for ORCA.
@@ -185,6 +374,8 @@ def _gen_input(
         orca_input += f"! DEFGRID{self.cfg.gridsize}\n"
         orca_input += "! MiniPrint\n"
         orca_input += "! NoTRAH\n"
+        if use_xtb_driver:
+            orca_input += "! Engrad\n"
         # "! AutoAux" keyword for super-heavy elements as def2/J ends at Rn
         if any(atom >= 86 for atom in molecule.ati):
             orca_input += "! AutoAux\n"