diff --git a/CHANGELOG.md b/CHANGELOG.md index c4f7308..805df4d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Distance constraints for xTB (normalized parsing, config validation, runtime enforcement, and optional integration test). +- Distance constrains for ORCA (uses xTB as as driver to set the same distance constrains for the postprocess). ## [0.6.0] - 2025-04-01 diff --git a/mindlessgen.toml b/mindlessgen.toml index 217cbe4..2d4ab68 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -106,6 +106,8 @@ basis = "def2-SVP" gridsize = 1 # > Maximum number of SCF cycles: Options: scf_cycles = 100 +# > Use xTB as an external ORCA driver to keep distance constraints during postprocessing. Options: +use_xtb_driver = false [turbomole] # > Path to the ridft executable. The name `ridft` is automatically searched for. Options: diff --git a/src/mindlessgen/generator/main.py b/src/mindlessgen/generator/main.py index 7087ce0..dcfc571 100644 --- a/src/mindlessgen/generator/main.py +++ b/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) diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index efb49fb..e01590d 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -1135,6 +1135,7 @@ def __init__(self: ORCAConfig) -> None: self._basis: str = "def2-SVP" self._gridsize: int = 1 self._scf_cycles: int = 100 + self._use_xtb_driver: bool = False def get_identifier(self) -> str: return "orca" @@ -1225,6 +1226,22 @@ def scf_cycles(self, max_scf_cycles: int): raise ValueError("Max SCF cycles should be greater than 0.") self._scf_cycles = max_scf_cycles + @property + def use_xtb_driver(self) -> bool: + """ + Determine whether the xTB external driver can be used during post-processing. + """ + return self._use_xtb_driver + + @use_xtb_driver.setter + def use_xtb_driver(self, enabled: bool): + """ + Enable or disable the usage of the xTB external driver during post-processing. + """ + if not isinstance(enabled, bool): + raise TypeError("use_xtb_driver should be a boolean.") + self._use_xtb_driver = enabled + class TURBOMOLEConfig(BaseConfig): """ diff --git a/src/mindlessgen/qm/orca.py b/src/mindlessgen/qm/orca.py index 0516ced..041bef0 100644 --- a/src/mindlessgen/qm/orca.py +++ b/src/mindlessgen/qm/orca.py @@ -8,8 +8,9 @@ from tempfile import TemporaryDirectory from ..molecules import Molecule -from ..prog import ORCAConfig +from ..prog import ORCAConfig, XTBConfig from .base import QMMethod +from .xtb import XTB, get_xtb_path class ORCA(QMMethod): @@ -17,7 +18,9 @@ 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 +31,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 +55,27 @@ 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) + if self.cfg.use_xtb_driver: + optimized_molecule = self.optimize_xtb_driver( + temp_path=temp_path, + molecule=molecule, + xyz_filename=xyz_filename, + ncores=ncores, + max_cycles=max_cycles, + verbosity=verbosity, + ) + return optimized_molecule inputname = "orca_opt.inp" orca_input = self._gen_input( - molecule, "molecule.xyz", ncores, True, max_cycles + molecule, + xyz_filename, + temp_path, + ncores, + True, + max_cycles, ) if verbosity > 1: print("ORCA input file:\n##################") @@ -63,12 +83,10 @@ def optimize( print("##################") with open(temp_path / inputname, "w", encoding="utf8") as f: f.write(orca_input) - # run orca arguments = [ inputname, ] - orca_log_out, orca_log_err, return_code = self._run( temp_path=temp_path, arguments=arguments ) @@ -78,7 +96,6 @@ def optimize( raise RuntimeError( f"ORCA failed with return code {return_code}:\n{orca_log_err}" ) - # read the optimized molecule from the output file xyzfile = Path(temp_path / inputname).resolve().with_suffix(".xyz") optimized_molecule = molecule.copy() @@ -103,10 +120,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) @@ -174,6 +191,7 @@ def _gen_input( self, molecule: Molecule, xyzfile: str, + _temp_path: Path, ncores: int, optimization: bool = False, opt_cycles: int | None = None, @@ -200,6 +218,145 @@ def _gen_input( orca_input += f"* xyzfile {molecule.charge} {molecule.uhf + 1} {xyzfile}\n" return orca_input + def optimize_xtb_driver( + self, + temp_path: Path, + molecule: Molecule, + xyz_filename: str, + ncores: int, + max_cycles: int | None = None, + verbosity: int = 1, + ) -> Molecule: + """ + Optimize a molecule using ORCA through the xTB external driver. + """ + + xtb_input = temp_path / "xtb.inp" + inputname = "orca_opt.inp" + self._write_xtb_input(molecule, xtb_input, inputname) + orca_input = self._gen_input_xtb_driver( + molecule, + xyz_filename, + temp_path, + ncores, + True, + max_cycles, + ) + if verbosity > 1: + print("ORCA input file:\n##################") + print(orca_input) + print("##################") + print("XTB input file:\n##################") + print(xtb_input) + print("##################") + with open(temp_path / inputname, "w", encoding="utf8") as f: + f.write(orca_input) + # run orca with xTB as a 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, + ) + if verbosity > 2: + print(orca_log_out) + if return_code != 0: + raise RuntimeError( + f"ORCA failed with return code {return_code}:\n{orca_log_err}" + ) + + # read the optimized molecule from the output file + xyzfile = temp_path / "xtbopt.xyz" + if not xyzfile.exists(): + raise RuntimeError( + "xTB-driven ORCA optimization did not produce 'xtbopt.xyz'." + ) + optimized_molecule = molecule.copy() + optimized_molecule.read_xyz_from_file(xyzfile) + return optimized_molecule + + 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 = get_xtb_path() + if self.xtb_cfg is None: + raise RuntimeError( + "xTB driver requested but no xTB configuration provided." + ) + xtb_runner = XTB(path=xtb_executable, xtb_config=self.xtb_cfg) + arguments = [ + 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]) + xtb_log_out, xtb_log_err, returncode = xtb_runner._run( + temp_path=temp_path, arguments=arguments + ) + return xtb_log_out, xtb_log_err, returncode + + 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." + ) + xtb_path = get_xtb_path() + xtb_writer = XTB(xtb_path, self.xtb_cfg) + generated = xtb_writer._prepare_distance_constraint_file( + molecule, xtb_input.parent + ) + if not generated: + raise RuntimeError( + "xTB driver requested but no distance constraints were generated." + ) + with xtb_input.open("a", encoding="utf8") as handle: + handle.write("$external\n") + handle.write(f" orca input file= {input_file}\n") + handle.write(f" orca bin= {self.path}\n") + handle.write("$end\n") + + def _gen_input_xtb_driver( + self, + molecule: Molecule, + xyzfile: str, + temp_path: Path, + ncores: int, + optimization: bool = False, + opt_cycles: int | None = None, + ) -> str: + """ + Generate a default input file for ORCA. + """ + orca_input = f"! {self.cfg.functional} {self.cfg.basis}\n" + orca_input += f"! DEFGRID{self.cfg.gridsize}\n" + orca_input += "! MiniPrint\n" + orca_input += "! NoTRAH\n" + 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" + orca_input += f"%scf\n\tMaxIter {self.cfg.scf_cycles}\n" + if not optimization: + orca_input += "\tConvergence Medium\n" + orca_input += "end\n" + orca_input += f"%pal nprocs {ncores} end\n\n" + orca_input += f"* xyzfile {molecule.charge} {molecule.uhf + 1} {xyzfile}\n" + return orca_input + # TODO: 1. Convert this to a @staticmethod of Class ORCA # 2. Rename to `get_method` or similar to enable an abstract interface diff --git a/test/test_qm/test_orca.py b/test/test_qm/test_orca.py new file mode 100644 index 0000000..4ea1a29 --- /dev/null +++ b/test/test_qm/test_orca.py @@ -0,0 +1,244 @@ +import subprocess as sp +from pathlib import Path +from types import SimpleNamespace +import numpy as np +import pytest +from mindlessgen.molecules import Molecule # type: ignore +from mindlessgen.prog import DistanceConstraint, XTBConfig # type: ignore +from mindlessgen.qm.orca import ORCA + + +class DummyORCAConfig(SimpleNamespace): + def __init__(self, **kwargs): + defaults = dict( + functional="B3LYP", + basis="def2-SVP", + gridsize=2, + scf_cycles=50, + optlevel="", + xtb_driver_path=None, + xtb_path=None, + use_xtb_driver=False, + ) + defaults.update(kwargs) + super().__init__(**defaults) + + +@pytest.fixture +def make_orca(): + def _factory(cfg=None, xtb_cfg_param=None): + cfg = cfg or DummyORCAConfig() + return ORCA(path="/usr/bin/orca", orcacfg=cfg, xtb_config=xtb_cfg_param) + + return _factory + + +@pytest.fixture +def xtb_cfg_oh(): + """ + xTB config with an O-H distance constraint. + """ + cfg = XTBConfig() + cfg.distance_constraints = [DistanceConstraint.from_cli_string("O,H,1.0")] + cfg.distance_constraint_force_constant = 0.5 + return cfg + + +@pytest.fixture +def xtb_cfg_ff(): + """ + xTB config with an F-F distance constraint. + """ + cfg = XTBConfig() + cfg.distance_constraints = [DistanceConstraint.from_cli_string("F,F,1.0")] + return cfg + + +@pytest.fixture +def xtb_cfg_fe3(): + """ + xTB config with a Fe-Fe distance constraint. + """ + cfg = XTBConfig() + cfg.distance_constraints = [ + DistanceConstraint.from_mapping({"pair": ["Fe", "Fe"], "distance": 2.5}) + ] + return cfg + + +@pytest.fixture +def mol_oh(): + """ + Simple O-H molecule for constraint tests. + """ + mol = Molecule("OH") + mol.ati = np.array([7, 0]) + return mol + + +@pytest.fixture +def mol_h2(): + """ + Simple H2 molecule for constraint tests. + """ + mol = Molecule("H2") + mol.ati = np.array([0, 0]) + return mol + + +@pytest.fixture +def mol_fe3(): + """ + Simple Fe3 molecule for constraint tests. + """ + mol = Molecule("Fe3") + mol.ati = np.array([25, 25, 25]) + return mol + + +@pytest.fixture +def fake_xtb_path(monkeypatch): + """ + Force xTB path discovery to a fake binary. + """ + xtb_path = Path("/fake/xtb") + monkeypatch.setattr("mindlessgen.qm.orca.get_xtb_path", lambda: xtb_path) + return xtb_path + + +def test_run_xtb_driver_success(monkeypatch, tmp_path, make_orca, fake_xtb_path): + orca = make_orca( + cfg=DummyORCAConfig(optlevel="tight"), + xtb_cfg_param=XTBConfig(), + ) + captured = {} + + def fake_run(args, cwd, capture_output, check): + captured["args"] = args + assert cwd == tmp_path + assert capture_output and check + return SimpleNamespace(stdout=b"ok", stderr=b"") + + monkeypatch.setattr("mindlessgen.qm.xtb.sp.run", fake_run) + out, err, code = orca._run_xtb_driver(tmp_path, "geom.xyz", "ctrl.inp", ncores=4) + assert captured["args"] == [ + str(fake_xtb_path), + "geom.xyz", + "--opt", + "tight", + "--orca", + "-I", + "ctrl.inp", + ] + assert out == "ok" + assert err == "" + assert code == 0 + + +def test_run_xtb_driver_failure_returns_error( + monkeypatch, tmp_path, make_orca, fake_xtb_path +): + """Ensure the ORCA wrapper surfaces errors from the xTB driver.""" + orca = make_orca(xtb_cfg_param=XTBConfig()) + + def fake_run(*_, **kwargs): + del kwargs + raise sp.CalledProcessError(1, "xtb", output=b"bad", stderr=b"worse") + + monkeypatch.setattr("mindlessgen.qm.xtb.sp.run", fake_run) + out, err, code = orca._run_xtb_driver( # pylint: disable=protected-access + tmp_path, "geom.xyz", "ctrl.inp", ncores=1 + ) + assert (out, err, code) == ("bad", "worse", 1) + + +def test_run_xtb_driver_requires_xtb_cfg( + monkeypatch, tmp_path, make_orca, fake_xtb_path +): + orca = make_orca() + with pytest.raises(RuntimeError, match="xTB driver requested"): + orca._run_xtb_driver(tmp_path, "geom.xyz", "ctrl.inp", ncores=1) + + +def test_write_xtb_input_creates_expected_file( + monkeypatch, tmp_path, make_orca, fake_xtb_path, mol_oh +): + xtb_cfg_instance = XTBConfig() + xtb_cfg_instance.distance_constraints = [] + xtb_cfg_instance.distance_constraint_force_constant = 0.7 + orca = make_orca(xtb_cfg_param=xtb_cfg_instance) + + def fake_prepare(self, molecule, temp_dir): + assert temp_dir == tmp_path + (temp_dir / "xtb.inp").write_text( + "\n".join( + [ + "$constrain", + " force constant= 0.7", + " distance: 1, 2, 1.00000", + "$end", + "", + ] + ), + encoding="utf8", + ) + return True + + monkeypatch.setattr( + "mindlessgen.qm.orca.XTB._prepare_distance_constraint_file", fake_prepare + ) + target = tmp_path / "xtb.inp" + orca._write_xtb_input(mol_oh, target, "orca.inp") + content = target.read_text().splitlines() + assert content[:4] == [ + "$constrain", + " force constant= 0.7", + " distance: 1, 2, 1.00000", + "$end", + ] + assert "$external" in content + assert " orca input file= orca.inp" in content + assert f" orca bin= {orca.path}" in content + + +def test_write_xtb_input_generates_constraints( + fake_xtb_path, tmp_path, make_orca, xtb_cfg_oh, mol_oh +): + orca = make_orca(xtb_cfg_param=xtb_cfg_oh) + + xtb_input = tmp_path / "xtb.inp" + orca._write_xtb_input(mol_oh, xtb_input, "orca.inp") + + contents = xtb_input.read_text(encoding="utf8").splitlines() + assert contents[0] == "$constrain" + assert "force constant= 0.5" in contents[1] + distance_lines = [line for line in contents if line.startswith(" distance:")] + assert distance_lines == [" distance: 1, 2, 1.0"] + assert "$external" in contents + assert " orca input file= orca.inp" in contents + assert f" orca bin= {orca.path}" in contents + + +def test_write_xtb_input_missing_atoms( + fake_xtb_path, tmp_path, make_orca, xtb_cfg_ff, mol_h2 +): + orca = make_orca(xtb_cfg_param=xtb_cfg_ff) + + with pytest.raises(RuntimeError): + orca._write_xtb_input(mol_h2, tmp_path / "xtb.inp", "orca.inp") + + +def test_distance_constraints_use_first_atoms( + fake_xtb_path, tmp_path, make_orca, xtb_cfg_fe3, mol_fe3 +): + orca = make_orca(xtb_cfg_param=xtb_cfg_fe3) + + xtb_input = tmp_path / "xtb.inp" + orca._write_xtb_input(mol_fe3, xtb_input, "orca.inp") + + contents = xtb_input.read_text(encoding="utf8").splitlines() + distance_lines = [line for line in contents if "distance:" in line] + + assert len(distance_lines) == 1 + assert "1, 2" in distance_lines[0] + assert ", 3," not in distance_lines[0]