diff --git a/Makefile b/Makefile index 08ebdd70eb..343134ac71 100644 --- a/Makefile +++ b/Makefile @@ -35,6 +35,7 @@ help: @echo " install-kinbot Install KinBot" @echo " install-sella Install Sella" @echo " install-xtb Install xTB" + @echo " install-crest Install CREST" @echo " install-torchani Install TorchANI" @echo " install-uma Install UMA (fairchem MLIP, gated model; users only, not CI)" @echo " install-ob Install OpenBabel" @@ -99,6 +100,9 @@ install-sella: install-xtb: bash $(DEVTOOLS_DIR)/install_xtb.sh +install-crest: + bash $(DEVTOOLS_DIR)/install_crest.sh + install-torchani: bash $(DEVTOOLS_DIR)/install_torchani.sh diff --git a/arc/constants.pxd b/arc/constants.pxd index 67b3f412bb..9fc3b9127d 100644 --- a/arc/constants.pxd +++ b/arc/constants.pxd @@ -1 +1 @@ -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 +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, angstrom_to_bohr diff --git a/arc/constants.py b/arc/constants.py index 771c09a05d..1386699e48 100644 --- a/arc/constants.py +++ b/arc/constants.py @@ -42,6 +42,7 @@ #: The Bohr radius :math:`a_0` in :math:`\mathrm{m}` a0 = 5.2917721092e-11 bohr_to_angstrom = a0 * 1e10 +angstrom_to_bohr = 1 / bohr_to_angstrom #: The atomic mass unit in :math:`\mathrm{kg}` amu = 1.660538921e-27 @@ -79,6 +80,7 @@ #: Vacuum permittivity epsilon_0 = 8.8541878128 + # 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 @@ -101,4 +103,5 @@ 'F': F, 'epsilon_0': epsilon_0, 'bohr_to_angstrom': bohr_to_angstrom, + 'angstrom_to_bohr': angstrom_to_bohr, }) diff --git a/arc/job/adapter.py b/arc/job/adapter.py index 4256b914e4..fa64d7bc8c 100644 --- a/arc/job/adapter.py +++ b/arc/job/adapter.py @@ -113,6 +113,7 @@ class JobEnum(str, Enum): autotst = 'autotst' # AutoTST, 10.1021/acs.jpca.7b07361, 10.26434/chemrxiv.13277870.v2 gcn = 'gcn' # Graph neural network for isomerization, https://doi.org/10.1021/acs.jpclett.0c00500 heuristics = 'heuristics' # ARC's heuristics + crest = 'crest' # CREST conformer/TS search kinbot = 'kinbot' # KinBot, 10.1016/j.cpc.2019.106947 linear = 'linear' # ARC's linear TS search goflow = 'goflow' # GoFlow, flow-matching E(3)-equivariant TS generator (Galustian et al., Digital Discovery 2025, 10.1039/D5DD00283D); https://github.com/heid-lab/goflow_lean diff --git a/arc/job/adapters/common.py b/arc/job/adapters/common.py index 43d567cddc..6271459e40 100644 --- a/arc/job/adapters/common.py +++ b/arc/job/adapters/common.py @@ -73,10 +73,10 @@ 'R_Addition_MultipleBond': ['autotst', 'kinbot', 'goflow', 'rits', 'linear'], 'Retroene': ['kinbot', 'goflow', 'rits', 'linear'], 'Singlet_Carbene_Intra_Disproportionation': ['gcn', 'xtb_gsm', 'orca_neb', 'goflow', 'rits', 'linear'], - 'XY_Addition_MultipleBond': ['goflow', 'rits', 'linear'], + 'XY_Addition_MultipleBond': ['heuristics', 'crest', 'goflow', 'rits', 'linear'], 'XY_elimination_hydroxyl': ['goflow', 'rits', 'linear'], - 'carbonyl_based_hydrolysis': ['heuristics'], - 'ether_hydrolysis': ['heuristics'], + 'carbonyl_based_hydrolysis': ['heuristics', 'crest'], + 'ether_hydrolysis': ['heuristics', 'crest'], 'halocarbene_recombination': ['goflow', 'rits', 'linear'], 'halocarbene_recombination_double': ['goflow', 'rits', 'linear'], 'intra_H_migration': ['autotst', 'gcn', 'kinbot', 'xtb_gsm', 'orca_neb', 'goflow', 'rits', 'linear'], @@ -88,7 +88,7 @@ 'intra_substitutionS_cyclization': ['goflow', 'rits', 'linear'], 'intra_substitutionS_isomerization': ['gcn', 'xtb_gsm', 'orca_neb', 'goflow', 'rits', 'linear'], 'lone_electron_pair_bond': ['goflow', 'rits', 'linear'], - 'nitrile_hydrolysis': ['heuristics'] + 'nitrile_hydrolysis': ['heuristics', 'crest'] } all_families_ts_adapters = [] diff --git a/arc/job/adapters/ts/__init__.py b/arc/job/adapters/ts/__init__.py index e1e5306438..31761c481e 100644 --- a/arc/job/adapters/ts/__init__.py +++ b/arc/job/adapters/ts/__init__.py @@ -1,4 +1,5 @@ import arc.job.adapters.ts.autotst_ts +import arc.job.adapters.ts.crest import arc.job.adapters.ts.gcn_ts import arc.job.adapters.ts.goflow_ts import arc.job.adapters.ts.heuristics @@ -6,4 +7,5 @@ import arc.job.adapters.ts.linear import arc.job.adapters.ts.orca_neb import arc.job.adapters.ts.rits_ts +import arc.job.adapters.ts.seed_hub import arc.job.adapters.ts.xtb_gsm diff --git a/arc/job/adapters/ts/crest.py b/arc/job/adapters/ts/crest.py new file mode 100644 index 0000000000..7554e0a5eb --- /dev/null +++ b/arc/job/adapters/ts/crest.py @@ -0,0 +1,881 @@ +""" +Utilities for running CREST within ARC. + +Separated from heuristics so CREST can be conditionally imported and reused. + +``TERMINAL_JOB_STATUSES`` lists every status after which a job is no longer polled. +:func:`arc.job.local.check_job_status` reports ``'done'``, ``'running'`` or ``'errored'`` and never +reports ``'failed'``, so ``'errored'`` has to be part of the terminal set for an errored job to stop +being polled. ``'failed'`` is retained because :func:`process_completed_jobs` reports it separately. + +``MAX_CREST_SEEDS`` bounds how many CREST jobs a single reaction spawns. Each seed becomes one +multi-core metadynamics job, and the adapter executes in core, so the seed list is capped to keep +both the number of concurrent jobs and the time ARC's main loop spends blocked on them bounded. +""" + +import datetime +import math +import os +import time + +import numpy as np +from typing import TYPE_CHECKING, List, Optional, Union + +from arc.common import almost_equal_coords, calc_rmsd, get_logger +from arc.exceptions import ConverterError +from arc.imports import settings, submit_scripts +from arc.job.adapter import JobAdapter +from arc.job.adapters.common import _initialize_adapter, ts_adapters_by_rmg_family +from arc.job.adapters.ts.heuristics import DIHEDRAL_INCREMENT +from arc.job.adapters.ts.seed_hub import get_backup_ts_seeds, get_ts_seeds, get_wrapper_constraints +from arc.job.factory import register_job_adapter +from arc.job.local import check_job_status, delete_job, submit_job +from arc.plotter import save_geo +from arc.species.converter import reorder_xyz_string, xyz_file_format_to_xyz, xyz_to_str +from arc.species.species import ARCSpecies, TSGuess, colliding_atoms + +if TYPE_CHECKING: + from arc.level import Level + from arc.reaction import ARCReaction + +logger = get_logger() + +MAX_CHECK_INTERVAL_SECONDS = 100 +TERMINAL_JOB_STATUSES = ('done', 'errored', 'failed') +DEFAULT_MAX_CREST_WALL_TIME = 24 * 60 * 60 +MAX_CREST_SEEDS = 8 + +# Number of atoms CREST holds rigid for each family it supports: the A--H--B triad of an +# H-abstraction, the *1/*2/*3/*4 ring of a four-centre XY addition, and the a/b/o/h1 tetrahedron of +# a hydrolysis. A family that is absent is never gated off by +# _crest_reactive_core_covers_molecule(). +CREST_REACTIVE_CORE_SIZES = {'H_Abstraction': 3, + 'XY_Addition_MultipleBond': 4, + 'carbonyl_based_hydrolysis': 4, + 'ether_hydrolysis': 4, + 'nitrile_hydrolysis': 4, + } + +CREST_PATH = settings.get("CREST_PATH", None) +CREST_ENV_PATH = settings.get("CREST_ENV_PATH", None) +SERVERS = settings.get("servers", {}) + + +def crest_available() -> bool: + """ + Return whether CREST is configured for use. + + CREST needs a configured local server and at least one of ``CREST_PATH`` (a standalone + executable) or ``CREST_ENV_PATH`` (an environment activation line). + + Returns: + bool: Whether CREST can be used. + """ + return bool(SERVERS.get("local")) and bool(CREST_PATH or CREST_ENV_PATH) + + +class CrestAdapter(JobAdapter): + """ + A class for executing CREST TS conformer searches based on heuristics-generated guesses. + """ + + def __init__(self, + project: str, + project_directory: str, + job_type: Union[List[str], str], + args: Optional[dict] = None, + bath_gas: Optional[str] = None, + checkfile: Optional[str] = None, + conformer: Optional[int] = None, + constraints: Optional[List] = None, + cpu_cores: Optional[str] = None, + dihedral_increment: Optional[float] = None, + dihedrals: Optional[List[float]] = None, + directed_scan_type: Optional[str] = None, + ess_settings: Optional[dict] = None, + ess_trsh_methods: Optional[List[str]] = None, + execution_type: Optional[str] = None, + fine: bool = False, + initial_time: Optional[Union['datetime.datetime', str]] = None, + irc_direction: Optional[str] = None, + job_id: Optional[int] = None, + job_memory_gb: float = 14.0, + job_name: Optional[str] = None, + job_num: Optional[int] = None, + job_server_name: Optional[str] = None, + job_status: Optional[List[Union[dict, str]]] = None, + level: Optional['Level'] = None, + max_job_time: Optional[float] = None, + run_multi_species: bool = False, + reactions: Optional[List['ARCReaction']] = None, + rotor_index: Optional[int] = None, + server: Optional[str] = None, + server_nodes: Optional[list] = None, + queue: Optional[str] = None, + attempted_queues: Optional[List[str]] = None, + species: Optional[List[ARCSpecies]] = None, + testing: bool = False, + times_rerun: int = 0, + torsions: Optional[List[List[int]]] = None, + tsg: Optional[int] = None, + xyz: Optional[dict] = None, + ): + + self.incore_capacity = 50 + self.job_adapter = 'crest' + self.command = None + self.execution_type = execution_type or 'incore' + + if reactions is None: + raise ValueError('Cannot execute TS CREST without ARCReaction object(s).') + + dihedral_increment = dihedral_increment or DIHEDRAL_INCREMENT + + _initialize_adapter(obj=self, + is_ts=True, + project=project, + project_directory=project_directory, + job_type=job_type, + args=args, + bath_gas=bath_gas, + checkfile=checkfile, + conformer=conformer, + constraints=constraints, + cpu_cores=cpu_cores, + dihedral_increment=dihedral_increment, + dihedrals=dihedrals, + directed_scan_type=directed_scan_type, + ess_settings=ess_settings, + ess_trsh_methods=ess_trsh_methods, + fine=fine, + initial_time=initial_time, + irc_direction=irc_direction, + job_id=job_id, + job_memory_gb=job_memory_gb, + job_name=job_name, + job_num=job_num, + job_server_name=job_server_name, + job_status=job_status, + level=level, + max_job_time=max_job_time, + run_multi_species=run_multi_species, + reactions=reactions, + rotor_index=rotor_index, + server=server, + server_nodes=server_nodes, + queue=queue, + attempted_queues=attempted_queues, + species=species, + testing=testing, + times_rerun=times_rerun, + torsions=torsions, + tsg=tsg, + xyz=xyz, + ) + + def write_input_file(self) -> None: + pass + + def set_files(self) -> None: + pass + + def set_additional_file_paths(self) -> None: + pass + + def set_input_file_memory(self) -> None: + pass + + def execute_incore(self): + """ + Run the CREST TS conformer search for every supported reaction of this job. + + For each reaction a seed list is built from the heuristics adapter, or, when that yields + nothing (e.g. a linear/cumulene reactive center such as HCCO that the heuristic Z-matrix + builder cannot assemble), from a successful TS guess that another adapter already appended + to ``rxn.ts_species.ts_guesses``. An external guess is a valid seed because CREST only needs + a seed geometry plus the family reactive-atom constraints, and those constraints are + re-derived from the seed geometry. Guesses whose method contains ``crest`` are excluded so + CREST is never seeded from a prior CREST result. + + The seed list is truncated to ``MAX_CREST_SEEDS`` entries, keeping the most geometrically + distinct ones, and the number of dropped seeds is logged. Each surviving seed is submitted + as one CREST job; the resulting geometries are appended as TS guesses unless they duplicate + an existing successful guess, in which case that guess records ``'crest'`` as an additional + method source. CREST is run without an explicit method flag, i.e., at its GFN2-xTB default. + """ + self._log_job_execution() + self.initial_time = self.initial_time if self.initial_time else datetime.datetime.now() + + supported_families = [key for key, val in ts_adapters_by_rmg_family.items() if 'crest' in val] + + self.reactions = [self.reactions] if not isinstance(self.reactions, list) else self.reactions + for rxn_index, rxn in enumerate(self.reactions): + if rxn.family not in supported_families: + logger.warning(f'The CREST TS search adapter does not support the {rxn.family} reaction family.') + continue + if any(spc.get_xyz() is None for spc in rxn.r_species + rxn.p_species): + logger.warning(f'The CREST TS search adapter cannot process a reaction if 3D coordinates of ' + f'some/all of its reactants/products are missing.\nNot processing {rxn}.') + continue + if not crest_available(): + logger.warning('CREST is not available. Skipping CREST TS search.') + break + + if _crest_reactive_core_covers_molecule(rxn): + logger.info( + f'Skipping CREST TS search for {rxn.label}: the reactive core spans essentially ' + f'the entire molecule (<=1 spectator atom), so CREST has no conformational degrees ' + f'of freedom to sample. Deferring to the other TS-search methods.' + ) + continue + + if rxn.ts_species is None: + rxn.ts_species = ARCSpecies(label='TS', + is_ts=True, + charge=rxn.charge, + multiplicity=rxn.multiplicity, + ) + + tsg = TSGuess(method='CREST') + tsg.tic() + + crest_job_dirs = [] + crest_references = {} + xyz_guesses = get_ts_seeds( + reaction=rxn, + base_adapter='heuristics', + dihedral_increment=self.dihedral_increment, + ) + if not xyz_guesses: + xyz_guesses = get_backup_ts_seeds(rxn, exclude_method='crest') + if xyz_guesses: + logger.info( + f'CREST heuristic seed construction failed for {rxn.label}; falling back to ' + f'{len(xyz_guesses)} external TS guess(es) as CREST seed(s).' + ) + if not xyz_guesses: + logger.warning(f'CREST TS search failed to generate any seed guesses for {rxn.label}.') + tsg.tok() + continue + + if len(xyz_guesses) > MAX_CREST_SEEDS: + logger.info( + f'The CREST seed generation produced {len(xyz_guesses)} seeds for {rxn.label}. ' + f'Submitting only the {MAX_CREST_SEEDS} most geometrically distinct of them and ' + f'dropping the other {len(xyz_guesses) - MAX_CREST_SEEDS}, so that a single reaction ' + f'cannot occupy the queue with an unbounded number of concurrent CREST jobs. Raise ' + f'MAX_CREST_SEEDS to sample more of them.' + ) + xyz_guesses = _select_diverse_seeds(seeds=xyz_guesses, max_seeds=MAX_CREST_SEEDS) + + for iteration, xyz_entry in enumerate(xyz_guesses): + xyz_guess = xyz_entry.get("xyz") + family = xyz_entry.get("family", rxn.family) + if xyz_guess is None: + continue + + crest_constraints = get_wrapper_constraints( + wrapper='crest', + reaction=rxn, + seed=xyz_entry, + ) + if not crest_constraints: + logger.warning( + f"Could not determine CREST constraint atoms for {rxn.label} crest seed {iteration} " + f"(family: {family}). Skipping this CREST seed." + ) + continue + + crest_job_dir = crest_ts_conformer_search( + xyz_guess, + constraints=crest_constraints, + path=self.local_path, + xyz_crest_int=f'{rxn_index}_{iteration}', + charge=rxn.charge if rxn.charge is not None else 0, + multiplicity=rxn.multiplicity if rxn.multiplicity is not None else 1, + ) + if crest_job_dir is None: + continue + crest_job_dirs.append(crest_job_dir) + crest_references[crest_job_dir] = { + 'xyz': xyz_guess, + 'constraints': crest_constraints, + } + + if not crest_job_dirs: + logger.warning(f'CREST TS search failed to prepare any jobs for {rxn.label}.') + tsg.tok() + continue + + crest_jobs = submit_crest_jobs(crest_job_dirs) + monitor_crest_jobs(crest_jobs) + xyz_guesses_crest = process_completed_jobs(crest_jobs, crest_references=crest_references) + tsg.tok() + + for method_index, xyz in enumerate(xyz_guesses_crest): + if xyz is None: + continue + unique = True + for other_tsg in rxn.ts_species.ts_guesses: + if not other_tsg.success or other_tsg.initial_xyz is None: + continue + if almost_equal_coords(xyz, other_tsg.initial_xyz): + if hasattr(other_tsg, "method_sources"): + other_tsg.method_sources = other_tsg._normalize_method_sources( + (other_tsg.method_sources or []) + ["crest"] + ) + unique = False + break + if unique: + ts_guess = TSGuess(method='CREST', + method_index=method_index, + t0=tsg.t0, + execution_time=tsg.execution_time, + success=True, + family=rxn.family, + xyz=xyz, + ) + rxn.ts_species.append_ts_guess(ts_guess) + save_geo(xyz=xyz, + path=self.local_path, + filename=f'CREST_{method_index}', + format_='xyz', + comment=f'CREST {method_index}, family: {rxn.family}', + ) + + if len(self.reactions) < 5: + successes = [tsg for tsg in rxn.ts_species.ts_guesses if tsg.success and 'crest' in tsg.method.lower()] + if successes: + logger.info(f'CREST successfully found {len(successes)} TS guesses for {rxn.label}.') + else: + logger.info(f'CREST did not find any successful TS guesses for {rxn.label}.') + + self.final_time = datetime.datetime.now() + + def execute_queue(self): + self.execute_incore() + + +def _crest_reactive_core_covers_molecule(rxn: 'ARCReaction') -> bool: + """ + Return whether the CREST reactive core spans essentially the whole molecule. + + The reactive core is the set of atoms CREST constrains: the three-center A--H--B triad of an + H_Abstraction, the four-center *1/*2/*3/*4 ring of an XY_Addition_MultipleBond, and the + four-center a/b/o/h1 tetrahedron of a hydrolysis. When the TS has at most one atom outside that + core there are no meaningful spectator conformational degrees of freedom for CREST metadynamics + to sample, so CREST cannot improve on the heuristic seed and should be skipped (the other + TS-search methods still cover the case). + + This is intentionally conservative: it only returns ``True`` for a family listed in + ``CREST_REACTIVE_CORE_SIZES`` that has <=1 spectator atom, never for a system that retains + real spectator degrees of freedom, and never for a family whose core size is unknown. + + Args: + rxn (ARCReaction): The reaction to check. + + Returns: + bool: Whether CREST should be skipped for this reaction. + """ + reactive_core_size = CREST_REACTIVE_CORE_SIZES.get(getattr(rxn, 'family', None)) + if reactive_core_size is None: + return False + reactant_species = getattr(rxn, 'r_species', None) or [] + atom_counts = [getattr(spc, 'number_of_atoms', None) for spc in reactant_species] + if not atom_counts or any(count is None for count in atom_counts): + return False + return (sum(atom_counts) - reactive_core_size) <= 1 + + +def _select_diverse_seeds(seeds: List[dict], max_seeds: int) -> List[dict]: + """ + Return at most ``max_seeds`` seed entries, preferring geometrically distinct ones. + + The first entry is always kept, and every further entry is the one whose smallest RMSD to the + already selected entries is the largest, so that near-duplicate seeds are dropped first. The + selected entries are returned in their original order. If the seed geometries cannot be compared + (a missing or differently shaped coordinate array), the first ``max_seeds`` entries are returned + instead. + + Args: + seeds (List[dict]): Seed entries in the schema returned by ``get_ts_seeds()``. + max_seeds (int): The maximum number of seed entries to return. + + Returns: + List[dict]: The selected seed entries. + """ + if max_seeds is None or max_seeds <= 0 or len(seeds) <= max_seeds: + return list(seeds) + coords = list() + for seed in seeds: + xyz = seed.get('xyz') if isinstance(seed, dict) else None + entry = None + if isinstance(xyz, dict) and xyz.get('coords'): + entry = np.array(xyz['coords'], dtype=float) + coords.append(entry) + if any(entry is None or entry.shape != coords[0].shape for entry in coords): + logger.debug('Could not compare the CREST seed geometries, keeping the first ' + f'{max_seeds} seeds.') + return list(seeds[:max_seeds]) + selected = [0] + remaining = list(range(1, len(seeds))) + while len(selected) < max_seeds and remaining: + best_index, best_distance = remaining[0], -1.0 + for index in remaining: + distance = min(calc_rmsd(coords[index], coords[chosen]) for chosen in selected) + if distance > best_distance: + best_index, best_distance = index, distance + selected.append(best_index) + remaining.remove(best_index) + return [seeds[index] for index in sorted(selected)] + + +def crest_ts_conformer_search( + xyz_guess: dict, + a_atom: Optional[int] = None, + h_atom: Optional[int] = None, + b_atom: Optional[int] = None, + path: str = "", + xyz_crest_int: Union[int, str] = 0, + constraints: Optional[dict] = None, + charge: int = 0, + multiplicity: int = 1, +) -> Optional[str]: + """ + Prepare a CREST TS conformer search job. + + Writes ``coords.ref`` (TURBOMOLE format, Bohr, ``x y z SYMBOL`` per line), ``constraints.inp`` + and a PBS/HTCondor submit script based on ``submit_scripts['local']['crest']`` into a + ``crest_`` sub-directory of ``path``, and returns that directory. A + ``crest_best.xyz`` left behind by an earlier run in the same directory is removed first, so a + run that dies on launch cannot have the previous run's geometry ingested as its own result. + If that file exists and cannot be removed, ``None`` is returned and no job is prepared: + :func:`process_completed_jobs` only tests whether ``crest_best.xyz`` exists, so leaving the + stale file in place would let it be reported as this run's result. + + ``constraints`` holds zero-based ``atoms`` and ``distance_pairs``, and, for the three-center + H-abstraction case, ``angle_atoms``. Every pair in ``distance_pairs`` is written as a + ``distance`` restraint; for the H-abstraction case the heavy--heavy A--B separation is written + as well. Those three distances rigidly determine the A--H--B triad, so no angle restraint is + emitted. The four-center hydrolysis case supplies all six pairwise distances among its reactive + atoms in ``distance_pairs`` and no ``angle_atoms``, which likewise determines its reactive core + without an angle restraint. Atoms outside the reactive zone form the ``$metadyn`` atom list; + when the reactive zone covers every atom no ``$metadyn`` block is written, since an atom-less + block only restates the CREST default. + + The restraint stiffness is written as ``force constant=0.5``: xtb parses ``key=value`` in this + block and silently ignores the colon form, falling back to its 0.05 default. + + ``--chrg`` and ``--uhf`` are appended to the CREST command for a charged or open-shell system; + without them CREST samples the neutral closed-shell species, which is the wrong electronic state + for the doublet TS of an H-abstraction. + + Args: + xyz_guess (dict): The seed geometry. + a_atom (int, optional): Legacy zero-based donor index, used only when ``constraints`` is None. + h_atom (int, optional): Legacy zero-based transferred-H index, used only when ``constraints`` is None. + b_atom (int, optional): Legacy zero-based acceptor index, used only when ``constraints`` is None. + path (str, optional): The parent directory in which the CREST job directory is created. + xyz_crest_int (Union[int, str], optional): The CREST job directory suffix, unique per seed. + constraints (dict, optional): The CREST constraint specification. + charge (int, optional): The net charge of the system. + multiplicity (int, optional): The spin multiplicity of the system. + + Returns: + Optional[str]: The CREST job directory path, or ``None`` if a stale ``crest_best.xyz`` + could not be removed or the local cluster software is not supported. + + Raises: + ValueError: If neither ``constraints`` nor a complete set of legacy A, H and B indices is given. + """ + if constraints is None: + if not all(isinstance(atom, int) for atom in (a_atom, h_atom, b_atom)): + raise ValueError('CREST requires either a constraint specification or legacy A, H, and B atom indices.') + constraints = { + 'atoms': (a_atom, h_atom, b_atom), + 'distance_pairs': ((a_atom, h_atom), (h_atom, b_atom)), + 'angle_atoms': (a_atom, h_atom, b_atom), + } + + path = os.path.join(path, f"crest_{xyz_crest_int}") + os.makedirs(path, exist_ok=True) + + stale_best_path = os.path.join(path, "crest_best.xyz") + if os.path.isfile(stale_best_path): + try: + os.remove(stale_best_path) + except OSError as e: + logger.error(f"Could not remove the stale CREST geometry at {stale_best_path}, skipping this CREST " + f"seed so that the previous run's geometry cannot be ingested as its own result: " + f"{e.__class__.__name__}: {e}") + return None + + symbols = xyz_guess["symbols"] + converted_coords = reorder_xyz_string( + xyz_str=xyz_to_str(xyz_guess), + reverse_atoms=True, + convert_to="bohr", + ) + coords_ref_content = f"$coord\n{converted_coords}\n$end\n" + coords_ref_path = os.path.join(path, "coords.ref") + with open(coords_ref_path, "w") as f: + f.write(coords_ref_content) + + num_atoms = len(symbols) + participating_atoms = tuple(atom + 1 for atom in constraints['atoms']) + distance_pairs = tuple((atom_1 + 1, atom_2 + 1) + for atom_1, atom_2 in constraints['distance_pairs']) + + angle_atoms = constraints.get('angle_atoms') + heavy_heavy_pair = None + if angle_atoms is not None: + heavy_heavy_pair = (angle_atoms[0] + 1, angle_atoms[2] + 1) + + list_of_atoms_numbers_not_participating_in_reaction = [ + i for i in range(1, num_atoms + 1) if i not in participating_atoms + ] + if not list_of_atoms_numbers_not_participating_in_reaction: + logger.warning(f"All {num_atoms} atoms are in the CREST reactive zone, so no $metadyn atom list " + f"is written and CREST will bias every atom by default.") + + constraints_path = os.path.join(path, "constraints.inp") + with open(constraints_path, "w") as f: + f.write("$constrain\n") + f.write(f" atoms: {', '.join(map(str, participating_atoms))}\n") + f.write(" force constant=0.5\n") + f.write(" reference=coords.ref\n") + for atom_1, atom_2 in distance_pairs: + f.write(f" distance: {atom_1}, {atom_2}, auto\n") + if heavy_heavy_pair is not None: + f.write(f" distance: {heavy_heavy_pair[0]}, {heavy_heavy_pair[1]}, auto\n") + if list_of_atoms_numbers_not_participating_in_reaction: + f.write("$metadyn\n") + f.write( + f' atoms: {", ".join(map(str, list_of_atoms_numbers_not_participating_in_reaction))}\n' + ) + f.write("$end\n") + + local_server = SERVERS.get("local", {}) + cpus = int(local_server.get("cpus", 8)) + if CREST_ENV_PATH: + crest_exe = "crest" + else: + crest_exe = CREST_PATH if CREST_PATH is not None else "crest" + + commands = [ + crest_exe, + "coords.ref", + "--cinp constraints.inp", + "--noreftopo", + f"-T {cpus}", + ] + if charge: + commands.append(f"--chrg {charge}") + if multiplicity and multiplicity > 1: + commands.append(f"--uhf {multiplicity - 1}") + command = " ".join(commands) + + activation_line = CREST_ENV_PATH or "" + + if SERVERS.get("local") is not None: + cluster_soft = SERVERS["local"]["cluster_soft"].lower() + local_templates = submit_scripts.get("local", {}) + crest_template = local_templates.get("crest") + crest_job_template = local_templates.get("crest_job") + + if cluster_soft in ["condor", "htcondor"]: + if crest_template is None: + crest_template = ( + "universe = vanilla\n" + "executable = job.sh\n" + "output = out.txt\n" + "error = err.txt\n" + "log = log.txt\n" + "request_cpus = {cpus}\n" + "request_memory = {memory}\n" + "JobBatchName = {name}\n" + "queue\n" + ) + if crest_job_template is None: + crest_job_template = ( + "#!/bin/bash -l\n" + "{activation_line}\n" + 'cd "{path}"\n' + "{commands}\n" + ) + sub_job = crest_template + format_params = { + "name": f"crest_{xyz_crest_int}", + "cpus": cpus, + "memory": int(SERVERS["local"].get("memory", 32.0) * 1024), + } + sub_job = sub_job.format(**format_params) + + with open( + os.path.join(path, settings["submit_filenames"]["HTCondor"]), "w" + ) as f: + f.write(sub_job) + + crest_job = crest_job_template.format( + path=path, + activation_line=activation_line, + commands=command, + ) + + with open(os.path.join(path, "job.sh"), "w") as f: + f.write(crest_job) + os.chmod(os.path.join(path, "job.sh"), 0o700) + + for fname in ("out.txt", "err.txt"): + fpath = os.path.join(path, fname) + if not os.path.exists(fpath): + with open(fpath, "w") as f: + f.write("") + os.chmod(fpath, 0o600) + + elif cluster_soft == "pbs": + if crest_template is None: + crest_template = ( + "#!/bin/bash -l\n" + "#PBS -q {queue}\n" + "#PBS -N {name}\n" + "#PBS -l select=1:ncpus={cpus}:mem={memory}gb\n" + "#PBS -o out.txt\n" + "#PBS -e err.txt\n\n" + "{activation_line}\n" + 'cd "{path}"\n' + "{commands}\n" + ) + sub_job = crest_template + format_params = { + "queue": SERVERS["local"].get("queue", "alon_q"), + "name": f"crest_{xyz_crest_int}", + "cpus": cpus, + "memory": int( + SERVERS["local"].get("memory", 32) + if SERVERS["local"].get("memory", 32) < 60 + else 40 + ), + "activation_line": activation_line, + "path": path, + "commands": command, + } + sub_job = sub_job.format(**format_params) + + submit_filename = settings["submit_filenames"]["PBS"] + submit_path = os.path.join(path, submit_filename) + with open(submit_path, "w") as f: + f.write(sub_job) + os.chmod(submit_path, 0o700) + + else: + logger.warning(f"CREST does not support the {cluster_soft!r} cluster software of the local server, " + f"skipping the CREST TS search.") + return None + + return path + + +def submit_crest_jobs(crest_paths: List[str]) -> dict: + """ + Submit CREST jobs to the server. + + A submission that yields no job id is skipped. ``submit_job()`` signals that either as + ``(None, None)`` or as ``('errored', '')``, the latter whenever the submit command wrote to + stderr or its stdout had an unrecognized shape, so the guard has to be on falsiness rather than + on ``None``. Keying the returned dict on such a value would collapse every failed submission + onto a single entry and lose the paths of all but the last one, and an empty-string key would in + addition adopt the first line of the queue listing as its status, because + ``check_job_status_in_stdout()`` matches a job id as a substring of the status line. + + Args: + crest_paths (List[str]): List of paths to the CREST directories. + + Returns: + dict: A dictionary keyed by job ID, holding the job path and status. + """ + crest_jobs = {} + for crest_path in crest_paths: + job_status, job_id = submit_job(path=crest_path) + if not job_id: + logger.error(f"Could not submit a CREST job for {crest_path} (got job ID {job_id!r} and status " + f"{job_status!r}), skipping it.") + continue + logger.debug(f"CREST job {job_id} submitted for {crest_path}") + crest_jobs[job_id] = {"path": crest_path, "status": job_status} + if crest_jobs: + job_ids = list(crest_jobs.keys()) + parent = os.path.dirname(crest_paths[0]) + logger.info(f"Submitted {len(job_ids)} CREST jobs ({job_ids[0]}-{job_ids[-1]}) for {parent}") + return crest_jobs + + +def monitor_crest_jobs(crest_jobs: dict, + check_interval: int = 300, + max_time_seconds: float = DEFAULT_MAX_CREST_WALL_TIME, + ) -> None: + """ + Monitor CREST jobs until they are complete. + + Polls every non-terminal job until all jobs reach a terminal status or ``max_time_seconds`` + elapses. A status query that raises leaves the job's status untouched and is retried on the next + poll, since a transient scheduler failure says nothing about the job. On expiry each + still-running job is cancelled with :func:`arc.job.local.delete_job` so + that it does not keep occupying the queue unowned; a cancellation failure is logged and does not + stop the remaining jobs from being cancelled. A cancelled job that had already written + ``crest_best.xyz`` is marked ``'done'`` so that geometry is still ingested downstream, and one + that wrote nothing is marked ``'errored'``. + + Args: + crest_jobs (dict): Dictionary containing job information (job ID, path, and status). + check_interval (int): Time interval (in seconds) to wait between status checks. + max_time_seconds (float): The wall time after which the remaining jobs are given up on. + """ + deadline = time.time() + max_time_seconds + while True: + all_done = True + for job_id, job_info in crest_jobs.items(): + if job_info["status"] not in TERMINAL_JOB_STATUSES: + try: + job_info["status"] = check_job_status(job_id) + except Exception as e: + logger.error(f"Could not check the status of the CREST job {job_id}, keeping its current " + f"status {job_info['status']!r} and retrying on the next poll: " + f"{e.__class__.__name__}: {e}") + if job_info["status"] not in TERMINAL_JOB_STATUSES: + all_done = False + if all_done: + break + if time.time() > deadline: + unfinished = {job_id: info["status"] for job_id, info in crest_jobs.items() + if info["status"] not in TERMINAL_JOB_STATUSES} + logger.error(f"CREST jobs did not finish within {max_time_seconds} s, cancelling them: {unfinished}.") + for job_id, info in crest_jobs.items(): + if info["status"] in TERMINAL_JOB_STATUSES: + continue + try: + delete_job(job_id) + except Exception as e: + logger.error(f"Could not cancel the timed-out CREST job {job_id}, it may still be occupying " + f"the queue: {e.__class__.__name__}: {e}") + crest_path = info.get("path") + crest_best_path = os.path.join(crest_path, "crest_best.xyz") if crest_path else None + if crest_best_path is not None and os.path.isfile(crest_best_path): + info["status"] = "done" + logger.info(f"The timed-out CREST job {job_id} had already written {crest_best_path}, " + f"ingesting that geometry.") + else: + info["status"] = "errored" + logger.info(f"The timed-out CREST job {job_id} wrote no geometry, dropping it.") + break + time.sleep(min(check_interval, MAX_CHECK_INTERVAL_SECONDS)) + + +def process_completed_jobs(crest_jobs: dict, crest_references: dict) -> list: + """ + Process the completed CREST jobs and update XYZ guesses. + + Only jobs with a ``'done'`` status are read. A geometry is rejected if it cannot be parsed + (CREST killed mid-write leaves a truncated file, which must not abort ARC), is empty, holds a + non-finite coordinate, has colliding atoms, has no reference seed recorded for its path, or no + longer preserves the seed's constrained reactive zone. + + Args: + crest_jobs (dict): Dictionary containing job information. + crest_references (dict): Reference seed geometry and constraint specification, keyed by CREST path. + + Returns: + list: The accepted CREST geometries. + """ + xyz_guesses = [] + for job_id, job_info in crest_jobs.items(): + crest_path = job_info["path"] + if job_info["status"] == "done": + crest_best_path = os.path.join(crest_path, "crest_best.xyz") + if os.path.exists(crest_best_path): + with open(crest_best_path, "r") as f: + content = f.read() + try: + xyz_guess = xyz_file_format_to_xyz(content) + except (ConverterError, IndexError, ValueError) as e: + logger.warning(f"Could not parse the CREST geometry at {crest_best_path}, skipping it: " + f"{e.__class__.__name__}: {e}") + continue + if xyz_guess is None or not xyz_guess.get('symbols'): + logger.warning(f"The CREST geometry at {crest_best_path} is empty, skipping it.") + continue + if any(not np.isfinite(coord) for xyz_row in xyz_guess['coords'] for coord in xyz_row): + logger.warning(f"The CREST geometry at {crest_best_path} contains non-finite coordinates, " + f"skipping it.") + continue + if colliding_atoms(xyz_guess): + logger.warning(f"Rejecting the CREST geometry from {crest_path}: it has colliding atoms.") + continue + reference = crest_references.get(crest_path) + if reference is None: + logger.warning(f"Rejecting unvalidated CREST geometry from {crest_path}: reference data is missing.") + continue + if not _preserves_reactive_constraints( + xyz=xyz_guess, + reference_xyz=reference['xyz'], + constraints=reference['constraints'], + ): + logger.warning( + f"Rejecting CREST geometry from {crest_path}: it does not preserve the reactive constraints." + ) + continue + xyz_guesses.append(xyz_guess) + else: + logger.error(f"crest_best.xyz not found in {crest_path}") + elif job_info["status"] == "failed": + logger.error(f"CREST job failed for {crest_path}") + + return xyz_guesses + + +def _preserves_reactive_constraints(xyz: dict, reference_xyz: dict, constraints: dict) -> bool: + """Return whether a CREST geometry preserves the seed's constrained reactive zone.""" + symbols = xyz.get('symbols') if isinstance(xyz, dict) else None + reference_symbols = reference_xyz.get('symbols') if isinstance(reference_xyz, dict) else None + if not symbols or symbols != reference_symbols or not isinstance(constraints, dict): + return False + try: + atoms = tuple(constraints['atoms']) + distance_pairs = tuple(tuple(pair) for pair in constraints['distance_pairs']) + if (not atoms or len(set(atoms)) != len(atoms) + or any(not isinstance(atom, int) or not 0 <= atom < len(symbols) for atom in atoms) + or any(len(pair) != 2 or pair[0] not in atoms or pair[1] not in atoms + for pair in distance_pairs)): + return False + for atom_1, atom_2 in distance_pairs: + distance = math.dist(xyz['coords'][atom_1], xyz['coords'][atom_2]) + reference_distance = math.dist(reference_xyz['coords'][atom_1], reference_xyz['coords'][atom_2]) + if abs(distance - reference_distance) > max(0.5, reference_distance * 0.5): + return False + angle_atoms = constraints.get('angle_atoms') + if angle_atoms is not None: + atom_1, vertex, atom_2 = angle_atoms + angle = _get_angle(xyz['coords'][atom_1], xyz['coords'][vertex], xyz['coords'][atom_2]) + reference_angle = _get_angle( + reference_xyz['coords'][atom_1], + reference_xyz['coords'][vertex], + reference_xyz['coords'][atom_2], + ) + if abs(angle - reference_angle) > 30.0: + return False + except (IndexError, KeyError, TypeError, ValueError): + return False + return True + + +def _get_angle(point_1: tuple, vertex: tuple, point_2: tuple) -> float: + """Return the angle in degrees for ``point_1``--``vertex``--``point_2``.""" + vector_1 = tuple(coord - center for coord, center in zip(point_1, vertex)) + vector_2 = tuple(coord - center for coord, center in zip(point_2, vertex)) + norm_product = math.sqrt(sum(coord ** 2 for coord in vector_1) * sum(coord ** 2 for coord in vector_2)) + if not norm_product: + raise ValueError('Cannot determine an angle from a zero-length vector.') + cosine = sum(coord_1 * coord_2 for coord_1, coord_2 in zip(vector_1, vector_2)) / norm_product + return math.degrees(math.acos(max(-1.0, min(1.0, cosine)))) + +register_job_adapter('crest', CrestAdapter) diff --git a/arc/job/adapters/ts/crest_test.py b/arc/job/adapters/ts/crest_test.py new file mode 100644 index 0000000000..8a26abd204 --- /dev/null +++ b/arc/job/adapters/ts/crest_test.py @@ -0,0 +1,1178 @@ +#!/usr/bin/env python3 +# encoding: utf-8 + +""" +Unit tests for arc.job.adapters.ts.crest +""" + +import os +import tempfile +import types +import unittest +from unittest.mock import patch + +from arc import constants +from arc.job.adapters.common import ts_adapters_by_rmg_family +from arc.job.adapters.ts import crest as crest_mod +from arc.job.adapters.ts.crest import CrestAdapter +from arc.job.adapters.ts.seed_hub import get_ts_seeds, get_wrapper_constraints +from arc.reaction import ARCReaction +from arc.settings.settings import ts_adapters +from arc.species.converter import str_to_xyz, xyz_to_str +from arc.species.species import ARCSpecies, TSGuess + +WATER_XYZ = str_to_xyz("""O 0.0 0.0 0.0 + H 0.0 0.0 0.96 + H 0.9 0.0 0.0""") + +OH_OH_XYZ = str_to_xyz("""O 0.00000000 -0.02752832 -1.20590500 + H 0.00000000 -0.02752832 -0.03383145 + O 0.00000000 -0.02752832 1.12142787 + H 0.00000000 0.90131726 1.37454478""") + +OH_OH_CONSTRAINTS = {'atoms': (0, 1, 2), + 'distance_pairs': ((0, 1), (1, 2)), + 'angle_atoms': (0, 1, 2), + } + +CH4_OH_SEED_XYZ = str_to_xyz("""C 0.00000000 0.00000000 0.00000000 + H 0.00000000 0.00000000 1.30000000 + H 1.03000000 0.00000000 -0.36000000 + H -0.51000000 -0.89000000 -0.36000000 + H -0.51000000 0.89000000 -0.36000000 + O 0.00000000 0.00000000 2.58000000 + H 0.90000000 0.00000000 2.90000000""") + +CH4_OH_RESULT_XYZ = str_to_xyz("""C 0.01000000 0.00000000 0.00000000 + H 0.00000000 0.01000000 1.31000000 + H 1.04000000 0.00000000 -0.36000000 + H -0.51000000 -0.90000000 -0.36000000 + H -0.51000000 0.90000000 -0.35000000 + O 0.00000000 0.01000000 2.59000000 + H 0.91000000 0.00000000 2.90000000""") + +CH4_OH_SEED = {'xyz': CH4_OH_SEED_XYZ, + 'family': 'H_Abstraction', + 'method': 'Heuristics', + 'source_adapter': 'heuristics', + 'metadata': {'reactive_atoms': {'A': 0, 'H': 1, 'B': 5}}, + } + + +HYDROLYSIS_FAMILIES = ('carbonyl_based_hydrolysis', 'ether_hydrolysis', 'nitrile_hydrolysis') + +# A methyl formate + water four-center hydrolysis TS seed, as heuristics.hydrolysis() generates it. +# Its atom roles are a=2 (the carbonyl carbon), b=1 (the leaving ester oxygen), e=3, o=8 (the water +# oxygen), h1=9 (the transferring water hydrogen) and d=7; atom 10 is the spectator water hydrogen. +HYDROLYSIS_SEED_XYZ = str_to_xyz("""C 0.14934638 -1.04077992 -1.29688828 + O 0.14934638 -1.04077992 0.13108654 + C 0.14934638 0.57218659 0.84211324 + O -0.41188584 1.56260964 0.40276657 + H -0.83434378 -0.75219743 -1.68052037 + H 0.36843817 -2.05708179 -1.63593228 + H 0.92676740 -0.37254154 -1.68052028 + H 0.14934645 0.47946490 1.93926498 + O 1.84045116 -0.06034836 1.24530570 + H 1.33683018 -0.95702439 0.60778079 + H 1.85273271 -0.25279783 2.19594358""") + +HYDROLYSIS_POSITIONAL_INDICES = [2, 1, 3, 8, 9, 7] + +HYDROLYSIS_REACTIVE_ATOMS = {'a': 2, 'b': 1, 'o': 8, 'h1': 9} + +# All six pairwise separations among a, b, o and h1, zero-based. +HYDROLYSIS_DISTANCE_PAIRS = ((2, 1), (2, 8), (2, 9), (1, 8), (1, 9), (8, 9)) + + +def make_ch4_oh_reaction() -> ARCReaction: + """Return a CH4 + OH <=> CH3 + H2O ARCReaction.""" + return ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2O', smiles='O')], + ) + + +def make_hydrolysis_seed(family: str = 'carbonyl_based_hydrolysis', + indices=None, + reactive_atoms: bool = True, + ) -> dict: + """ + Return a hydrolysis seed entry in the schema that get_ts_seeds() produces. + + Args: + family (str, optional): The hydrolysis family to label the seed with. + indices (optional): The raw ``indices`` metadata entry. Defaults to the positional + ``(a, b, e, o, h1, d)`` sequence that hydrolysis() emits. + reactive_atoms (bool, optional): Whether to also carry the named ``reactive_atoms`` dict. + + Returns: + dict: The seed entry. + """ + metadata = {'indices': list(HYDROLYSIS_POSITIONAL_INDICES) if indices is None else indices} + if reactive_atoms: + metadata['reactive_atoms'] = dict(HYDROLYSIS_REACTIVE_ATOMS) + return {'xyz': HYDROLYSIS_SEED_XYZ, + 'family': family, + 'method': 'Heuristics', + 'source_adapter': 'heuristics', + 'metadata': metadata, + } + + +def make_methyl_formate_hydrolysis_reaction() -> ARCReaction: + """Return a methyl formate + water <=> formic acid + methanol ARCReaction.""" + return ARCReaction(r_species=[ARCSpecies(label='ester', smiles='COC=O'), + ARCSpecies(label='H2O', smiles='O')], + p_species=[ARCSpecies(label='formicacid', smiles='C(=O)O'), + ARCSpecies(label='methanol', smiles='CO')], + family='carbonyl_based_hydrolysis', + ) + + +class CrestTestCase(unittest.TestCase): + """A base class providing a temporary directory and CREST module-global patching.""" + + def setUp(self): + """Create a temporary directory that is removed after the test.""" + self.tmpdir = tempfile.TemporaryDirectory() + self.addCleanup(self.tmpdir.cleanup) + + def patch_crest_module(self, **overrides): + """ + Patch the CREST module globals for the duration of the test. + + Args: + overrides: Module global names mapped to the values to patch them with. Any global that + is not overridden gets a local-PBS default suitable for input generation. + + Returns: + module: The patched ``arc.job.adapters.ts.crest`` module. + """ + values = {'settings': {'submit_filenames': {'PBS': 'submit.sh', 'HTCondor': 'submit.sub'}}, + 'submit_scripts': {'local': {}}, + 'CREST_PATH': '/usr/bin/crest', + 'CREST_ENV_PATH': '', + 'SERVERS': {'local': {'cluster_soft': 'pbs', 'cpus': 4, 'memory': 8, 'queue': 'testq'}}, + } + values.update(overrides) + for name, value in values.items(): + patcher = patch.object(crest_mod, name, value) + patcher.start() + self.addCleanup(patcher.stop) + return crest_mod + + +class TestCrestInputGeneration(CrestTestCase): + """Tests for CREST input generation.""" + + def test_creates_valid_input_files(self): + """Ensure CREST inputs are written with the expected content and format.""" + self.patch_crest_module(submit_scripts={'local': { + 'crest': ('#PBS -q {queue}\n' + '#PBS -N {name}\n' + '#PBS -l select=1:ncpus={cpus}:mem={memory}gb\n'), + 'crest_job': '{activation_line}\ncd {path}\n{commands}\n', + }}) + + crest_dir = crest_mod.crest_ts_conformer_search(WATER_XYZ, 0, 1, 2, self.tmpdir.name, 0) + + coords_path = os.path.join(crest_dir, 'coords.ref') + constraints_path = os.path.join(crest_dir, 'constraints.inp') + submit_path = os.path.join(crest_dir, 'submit.sh') + + self.assertTrue(os.path.exists(coords_path)) + self.assertTrue(os.path.exists(constraints_path)) + self.assertTrue(os.path.exists(submit_path)) + + with open(coords_path) as f: + coords = f.read().strip().splitlines() + self.assertEqual(coords[0].strip(), '$coord') + self.assertEqual(coords[-1].strip(), '$end') + self.assertEqual(len(coords) - 2, len(WATER_XYZ['symbols'])) + + with open(constraints_path) as f: + constraints = f.read() + self.assertIn('atoms: 1, 2, 3', constraints) + self.assertIn('force constant=0.5', constraints) + self.assertIn('reference=coords.ref', constraints) + self.assertIn('distance: 1, 2, auto', constraints) + self.assertIn('distance: 2, 3, auto', constraints) + self.assertTrue(constraints.strip().endswith('$end')) + + def test_coords_ref_is_in_bohr_with_turbomole_column_order(self): + """coords.ref holds ``x y z SYMBOL`` lines whose coordinates are in Bohr, not Angstrom.""" + self.patch_crest_module() + + crest_dir = crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + constraints={'atoms': (0, 1, 2), + 'distance_pairs': ((0, 1), (1, 2))}, + path=self.tmpdir.name, + xyz_crest_int=0, + ) + with open(os.path.join(crest_dir, 'coords.ref')) as f: + lines = f.read().strip().splitlines() + atom_lines = lines[1:-1] + self.assertEqual(len(atom_lines), 3) + + parsed = list() + for line, symbol in zip(atom_lines, WATER_XYZ['symbols']): + tokens = line.split() + self.assertEqual(len(tokens), 4) + self.assertEqual(tokens[3], symbol) + parsed.append(tuple(float(token) for token in tokens[:3])) + + for parsed_coords, expected_coords in zip(parsed, WATER_XYZ['coords']): + for parsed_coord, expected_coord in zip(parsed_coords, expected_coords): + self.assertAlmostEqual(parsed_coord, expected_coord * constants.angstrom_to_bohr, places=6) + + self.assertAlmostEqual(parsed[1][2], 0.96 * constants.angstrom_to_bohr, places=6) + self.assertNotAlmostEqual(parsed[1][2], 0.96, places=3) + + def test_charge_and_multiplicity_flags(self): + """--chrg and --uhf are written for a charged open-shell system and omitted otherwise.""" + self.patch_crest_module() + constraints = {'atoms': (0, 1, 2), 'distance_pairs': ((0, 1), (1, 2))} + + anion_dir = crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + constraints=constraints, + path=self.tmpdir.name, + xyz_crest_int=0, + charge=-1, + multiplicity=2, + ) + with open(os.path.join(anion_dir, 'submit.sh')) as f: + anion_submit = f.read() + self.assertIn('--chrg -1', anion_submit) + self.assertIn('--uhf 1', anion_submit) + + neutral_dir = crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + constraints=constraints, + path=self.tmpdir.name, + xyz_crest_int=1, + charge=0, + multiplicity=1, + ) + with open(os.path.join(neutral_dir, 'submit.sh')) as f: + neutral_submit = f.read() + self.assertNotIn('--chrg', neutral_submit) + self.assertNotIn('--uhf', neutral_submit) + + def test_h_abstraction_pins_three_distances_and_writes_no_angle(self): + """The H-abstraction $constrain block pins A--H, H--B and A--B, which fixes the triad.""" + self.patch_crest_module() + + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=OH_OH_XYZ, + constraints=OH_OH_CONSTRAINTS, + path=self.tmpdir.name, + xyz_crest_int=3, + ) + with open(os.path.join(crest_path, 'constraints.inp')) as f: + constraints_text = f.read() + self.assertIn('distance: 1, 2, auto', constraints_text) + self.assertIn('distance: 2, 3, auto', constraints_text) + self.assertIn('distance: 1, 3, auto', constraints_text) + self.assertNotIn('angle:', constraints_text) + self.assertIn('$metadyn\n atoms: 4\n', constraints_text) + + def test_no_metadyn_block_when_every_atom_is_reactive(self): + """An atom-less $metadyn block is not written; the reactive zone covering all atoms is logged.""" + self.patch_crest_module() + + with self.assertLogs('arc', level='WARNING') as log: + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + constraints={'atoms': (0, 1, 2), + 'distance_pairs': ((0, 1), (1, 2))}, + path=self.tmpdir.name, + xyz_crest_int=0, + ) + with open(os.path.join(crest_path, 'constraints.inp')) as f: + constraints_text = f.read() + self.assertNotIn('$metadyn', constraints_text) + self.assertIn('$constrain', constraints_text) + self.assertTrue(constraints_text.strip().endswith('$end')) + self.assertTrue(any('$metadyn' in message for message in log.output)) + + def test_stale_crest_best_is_removed_before_launching(self): + """A crest_best.xyz left by an earlier run must not survive into a new run's directory.""" + self.patch_crest_module() + stale_dir = os.path.join(self.tmpdir.name, 'crest_0') + os.makedirs(stale_dir) + stale_path = os.path.join(stale_dir, 'crest_best.xyz') + with open(stale_path, 'w') as f: + f.write(f"4\nstale geometry\n{xyz_to_str(OH_OH_XYZ)}\n") + + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=OH_OH_XYZ, + constraints=OH_OH_CONSTRAINTS, + path=self.tmpdir.name, + xyz_crest_int=0, + ) + self.assertEqual(crest_path, stale_dir) + self.assertFalse(os.path.exists(stale_path)) + + def test_a_stale_crest_best_that_cannot_be_removed_aborts_the_seed(self): + """A stale crest_best.xyz that survives must abort the seed rather than be run beside. + + process_completed_jobs() accepts a geometry on the strength of crest_best.xyz existing, so a + job prepared next to a surviving stale file would report the previous run's geometry as its + own result. + """ + self.patch_crest_module() + stale_dir = os.path.join(self.tmpdir.name, 'crest_0') + os.makedirs(stale_dir) + stale_path = os.path.join(stale_dir, 'crest_best.xyz') + with open(stale_path, 'w') as f: + f.write(f"4\nstale geometry\n{xyz_to_str(OH_OH_XYZ)}\n") + + with patch('os.remove', side_effect=OSError(13, 'Permission denied')), \ + self.assertLogs('arc', level='ERROR') as log: + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=OH_OH_XYZ, + constraints=OH_OH_CONSTRAINTS, + path=self.tmpdir.name, + xyz_crest_int=0, + ) + self.assertIsNone(crest_path) + self.assertTrue(os.path.exists(stale_path)) + self.assertFalse(os.path.exists(os.path.join(stale_dir, 'coords.ref'))) + self.assertFalse(os.path.exists(os.path.join(stale_dir, 'constraints.inp'))) + self.assertTrue(any('Could not remove the stale CREST geometry' in message for message in log.output)) + + def test_creates_submit_file_without_crest_templates(self): + """Fallback submit template generation works when submit.py has no CREST templates.""" + self.patch_crest_module() + + crest_dir = crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + constraints=OH_OH_CONSTRAINTS, + path=self.tmpdir.name, + xyz_crest_int=1, + ) + submit_path = os.path.join(crest_dir, 'submit.sh') + self.assertTrue(os.path.exists(submit_path)) + with open(submit_path) as f: + submit_text = f.read() + self.assertIn('#PBS -q testq', submit_text) + self.assertIn('coords.ref --cinp constraints.inp --noreftopo -T 4', submit_text) + + def test_creates_xy_distance_constraints_and_validates_completed_geometry(self): + """Write all three XY recipe distances and reject a geometry that loses one.""" + reference_xyz = str_to_xyz("""C 0.0000 0.0000 0.6670 + C 0.0000 0.0000 -0.6670 + H 0.0000 0.9210 1.2320 + H 0.0000 -0.9210 1.2320 + H 0.0000 0.9210 -1.2320 + H 0.0000 -0.9210 -1.2320 + Cl 0.0000 2.1000 -0.6670 + H 0.0000 1.6000 0.6670""") + constraints = {'atoms': (1, 0, 7, 6), + 'distance_pairs': ((1, 7), (0, 6), (7, 6)), + } + self.patch_crest_module() + + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=reference_xyz, + constraints=constraints, + path=self.tmpdir.name, + xyz_crest_int=2, + ) + with open(os.path.join(crest_path, 'constraints.inp')) as f: + constraints_text = f.read() + self.assertIn('atoms: 2, 1, 8, 7', constraints_text) + self.assertIn('distance: 2, 8, auto', constraints_text) + self.assertIn('distance: 1, 7, auto', constraints_text) + self.assertIn('distance: 8, 7, auto', constraints_text) + self.assertIn('atoms: 3, 4, 5, 6', constraints_text) + + crest_best_path = os.path.join(crest_path, 'crest_best.xyz') + jobs = {'123': {'path': crest_path, 'status': 'done'}} + references = {crest_path: {'xyz': reference_xyz, 'constraints': constraints}} + with open(crest_best_path, 'w') as f: + f.write(f"8\nCREST geometry\n{xyz_to_str(reference_xyz)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), [reference_xyz]) + + dissociated_xyz = dict(reference_xyz) + dissociated_coords = list(reference_xyz['coords']) + dissociated_coords[6] = (0.0, 8.0, -0.6670) + dissociated_xyz['coords'] = tuple(dissociated_coords) + with open(crest_best_path, 'w') as f: + f.write(f"8\nCREST geometry\n{xyz_to_str(dissociated_xyz)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), []) + + def test_unsupported_cluster_software_skips_instead_of_raising(self): + """An unsupported local queue must skip CREST, not abort the ARC run. + + crest_ts_conformer_search() is reached from job.execute(), which Scheduler.run_job calls + with no try/except, so raising there kills the whole project. + """ + self.patch_crest_module(SERVERS={'local': {'cluster_soft': 'Slurm', 'cpus': 4}}) + result = crest_mod.crest_ts_conformer_search(WATER_XYZ, + path=self.tmpdir.name, + constraints={'atoms': [0, 1], + 'distance_pairs': [(0, 1)]}, + ) + self.assertIsNone(result) + + def test_constructor_and_legacy_argument_validation(self): + """The adapter rejects a missing reaction, and the legacy call rejects incomplete indices.""" + with self.assertRaises(ValueError): + CrestAdapter(job_type='tsg', + reactions=None, + testing=True, + project='test_CrestAdapter', + project_directory=self.tmpdir.name, + ) + with self.assertRaises(ValueError): + crest_mod.crest_ts_conformer_search(xyz_guess=WATER_XYZ, + a_atom=None, + h_atom=1, + b_atom=2, + path=self.tmpdir.name, + ) + + +class TestCrestBackupSeeds(CrestTestCase): + """Tests for seeding CREST from TS guesses that other adapters produced.""" + + def test_get_backup_ts_seeds_selects_successful_non_crest_guesses(self): + """The backup seed picker skips CREST/failed guesses, prefers opt_xyz, and dedups.""" + from arc.job.adapters.ts.seed_hub import get_backup_ts_seeds + + seed_geom = OH_OH_XYZ + other_geom = str_to_xyz("""O 0.10000000 -0.02752832 -1.20590500 + H 0.00000000 -0.02752832 -0.03383145 + O 0.00000000 -0.02752832 1.12142787 + H 0.00000000 0.90131726 1.37454478""") + + def make_tsg(method, success, initial_xyz, opt_xyz=None): + return types.SimpleNamespace(method=method, success=success, + initial_xyz=initial_xyz, opt_xyz=opt_xyz) + + ts_guesses = [ + make_tsg('crest', True, other_geom), + make_tsg('autotst', False, other_geom), + make_tsg('autotst', True, other_geom, opt_xyz=seed_geom), + make_tsg('heuristics', True, seed_geom), + ] + reaction = types.SimpleNamespace(family='H_Abstraction', + ts_species=types.SimpleNamespace(ts_guesses=ts_guesses), + ) + + seeds = get_backup_ts_seeds(reaction, exclude_method='crest') + self.assertEqual(len(seeds), 1) + seed = seeds[0] + self.assertEqual(seed['xyz'], seed_geom) + self.assertEqual(seed['source_adapter'], 'autotst') + self.assertEqual(seed['family'], 'H_Abstraction') + self.assertEqual(seed['metadata'], {}) + + self.assertEqual(get_backup_ts_seeds(types.SimpleNamespace(family='H_Abstraction')), []) + empty_rxn = types.SimpleNamespace(family='H_Abstraction', + ts_species=types.SimpleNamespace(ts_guesses=[])) + self.assertEqual(get_backup_ts_seeds(empty_rxn), []) + + def test_backup_seed_drives_constraint_derivation_and_input_generation(self): + """An external TS guess seeds CREST: constraints are re-derived from its geometry.""" + from arc.job.adapters.ts.seed_hub import get_backup_ts_seeds, get_wrapper_constraints + + reaction = types.SimpleNamespace( + family='H_Abstraction', + ts_species=types.SimpleNamespace(ts_guesses=[ + types.SimpleNamespace(method='autotst', success=True, + initial_xyz=OH_OH_XYZ, opt_xyz=None), + ]), + ) + + seeds = get_backup_ts_seeds(reaction, exclude_method='crest') + self.assertEqual(len(seeds), 1) + constraints = get_wrapper_constraints(wrapper='crest', reaction=reaction, seed=seeds[0]) + self.assertIsNotNone(constraints) + self.assertIn('atoms', constraints) + self.assertIn('distance_pairs', constraints) + self.assertIn('angle_atoms', constraints) + + self.patch_crest_module() + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=seeds[0]['xyz'], + constraints=constraints, + path=self.tmpdir.name, + xyz_crest_int=0, + ) + self.assertTrue(os.path.exists(os.path.join(crest_path, 'coords.ref'))) + with open(os.path.join(crest_path, 'constraints.inp')) as f: + constraints_text = f.read() + self.assertIn('atoms:', constraints_text) + self.assertIn('reference=coords.ref', constraints_text) + self.assertIn('$metadyn', constraints_text) + + +class TestCrestAvailability(CrestTestCase): + """Tests for crest_available().""" + + def test_crest_available_matrix(self): + """CREST is available only with a local server and at least one CREST path configured.""" + with patch.object(crest_mod, 'SERVERS', {}), \ + patch.object(crest_mod, 'CREST_PATH', '/usr/bin/crest'), \ + patch.object(crest_mod, 'CREST_ENV_PATH', ''): + self.assertFalse(crest_mod.crest_available()) + + with patch.object(crest_mod, 'SERVERS', {'local': {'cluster_soft': 'pbs'}}), \ + patch.object(crest_mod, 'CREST_PATH', None), \ + patch.object(crest_mod, 'CREST_ENV_PATH', None): + self.assertFalse(crest_mod.crest_available()) + + with patch.object(crest_mod, 'SERVERS', {'local': {'cluster_soft': 'pbs'}}), \ + patch.object(crest_mod, 'CREST_PATH', '/usr/bin/crest'), \ + patch.object(crest_mod, 'CREST_ENV_PATH', None): + self.assertTrue(crest_mod.crest_available()) + + with patch.object(crest_mod, 'SERVERS', {'local': {'cluster_soft': 'pbs'}}), \ + patch.object(crest_mod, 'CREST_PATH', None), \ + patch.object(crest_mod, 'CREST_ENV_PATH', 'conda activate crest_env'): + self.assertTrue(crest_mod.crest_available()) + + +class TestCrestGating(CrestTestCase): + """Tests for the reactive-core gate and the seed cap.""" + + @staticmethod + def make_rxn(family, reactant_atom_counts): + """Return a minimal reaction stand-in carrying a family and reactant atom counts.""" + return types.SimpleNamespace( + family=family, + r_species=[types.SimpleNamespace(number_of_atoms=n) for n in reactant_atom_counts], + ) + + def test_tiny_system_gate_skips_crest_only_for_whole_molecule_core(self): + """CREST is gated off only when at most one atom sits outside the constrained reactive core.""" + self.assertTrue(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('H_Abstraction', [2, 2]))) + self.assertTrue(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('H_Abstraction', [1, 2]))) + # 5 atoms, i.e. 2 spectators, is the first case that must NOT be gated off. + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('H_Abstraction', [3, 2]))) + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('H_Abstraction', [5, 2]))) + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('H_Abstraction', [None, 2]))) + # The XY-addition core is the four *1-*4 atoms CREST constrains, so a 4- or 5-atom system + # has no spectator to sample, while C2H4 + HCl (8 atoms) has four and must not be gated off. + self.assertTrue( + crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('XY_Addition_MultipleBond', [2, 2])) + ) + self.assertTrue( + crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('XY_Addition_MultipleBond', [3, 2])) + ) + self.assertFalse( + crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('XY_Addition_MultipleBond', [6, 2])) + ) + self.assertFalse( + crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('XY_Addition_MultipleBond', [None, 2])) + ) + + def test_tiny_system_gate_covers_the_four_center_hydrolysis_core(self): + """A hydrolysis core of four atoms is gated off only when at most one spectator remains.""" + for family in HYDROLYSIS_FAMILIES: + # 5 atoms, i.e. 1 spectator outside the a/b/o/h1 core, is gated off. + self.assertTrue(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn(family, [2, 3])), + msg=f'{family} was not gated off at 5 atoms.') + self.assertTrue(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn(family, [1, 3]))) + # 6 atoms, e.g. HCN + H2O, is the smallest real substrate and must NOT be gated off. + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn(family, [3, 3])), + msg=f'{family} was wrongly gated off at 6 atoms.') + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn(family, [8, 3]))) + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn(family, [None, 3]))) + self.assertFalse(crest_mod._crest_reactive_core_covers_molecule(self.make_rxn('intra_H_migration', [2, 3]))) + + def test_select_diverse_seeds_keeps_distinct_geometries(self): + """The seed selector keeps geometrically distinct seeds and drops near-duplicates.""" + def seed(offset): + coords = tuple(tuple(coord + offset for coord in atom_coords) for atom_coords in WATER_XYZ['coords']) + return {'xyz': {'symbols': WATER_XYZ['symbols'], 'isotopes': WATER_XYZ['isotopes'], 'coords': coords}} + + seeds = [seed(0.0), seed(0.001), seed(0.002), seed(5.0)] + selected = crest_mod._select_diverse_seeds(seeds=seeds, max_seeds=2) + self.assertEqual(len(selected), 2) + self.assertIn(seeds[0], selected) + self.assertIn(seeds[3], selected) + + self.assertEqual(crest_mod._select_diverse_seeds(seeds=seeds, max_seeds=10), seeds) + + incomparable = [{'xyz': WATER_XYZ}, {'xyz': OH_OH_XYZ}, {'xyz': WATER_XYZ}] + self.assertEqual(crest_mod._select_diverse_seeds(seeds=incomparable, max_seeds=2), incomparable[:2]) + + +class TestCrestJobLifecycle(CrestTestCase): + """Test the parts of the CREST lifecycle that can hang or abort the whole ARC run.""" + + def test_monitor_terminates_on_an_errored_job(self): + """A queue reporting 'errored' must end the wait. + + check_job_status() reports 'done', 'running' or 'errored' and never 'failed', so a terminal + set of ('done', 'failed') left an errored job polled forever, blocking the scheduler. + """ + crest_jobs = {'1': {'path': os.path.join(self.tmpdir.name, 'crest_0'), 'status': 'running'}} + calls = {'n': 0} + + def fake_status(job_id): + calls['n'] += 1 + return 'errored' + + with patch.object(crest_mod, 'check_job_status', fake_status): + crest_mod.monitor_crest_jobs(crest_jobs, check_interval=0, max_time_seconds=1) + self.assertEqual(crest_jobs['1']['status'], 'errored') + # Exactly one poll: 'errored' is terminal. If it were not, the loop would spin until the + # wall-clock deadline and poll many times, which is the hang this guards against. + self.assertEqual(calls['n'], 1) + + def test_a_transient_status_check_failure_leaves_the_status_unchanged(self): + """A scheduler error while polling says nothing about the job, so its status must survive. + + Marking the job 'errored' would make it terminal, so it would stop being polled and its + geometry would never be ingested even though the job went on to finish normally. + """ + crest_jobs = {'1': {'path': os.path.join(self.tmpdir.name, 'crest_0'), 'status': 'running'}} + calls = {'n': 0} + + def fake_status(job_id): + calls['n'] += 1 + if calls['n'] == 1: + raise RuntimeError('qstat is temporarily unavailable') + return 'done' + + with patch.object(crest_mod, 'check_job_status', fake_status), \ + self.assertLogs('arc', level='ERROR') as log: + crest_mod.monitor_crest_jobs(crest_jobs, check_interval=0, max_time_seconds=60) + + self.assertEqual(crest_jobs['1']['status'], 'done') + self.assertEqual(calls['n'], 2) + self.assertTrue(any('Could not check the status of the CREST job 1' in message + for message in log.output)) + + def test_monitor_cancels_and_ingests_on_the_wall_time(self): + """A timed-out job is cancelled; a geometry it already wrote is still ingested.""" + finished_path = os.path.join(self.tmpdir.name, 'crest_0') + empty_path = os.path.join(self.tmpdir.name, 'crest_1') + os.makedirs(finished_path) + os.makedirs(empty_path) + with open(os.path.join(finished_path, 'crest_best.xyz'), 'w') as f: + f.write(f"4\nCREST geometry\n{xyz_to_str(OH_OH_XYZ)}\n") + crest_jobs = {'1': {'path': finished_path, 'status': 'running'}, + '2': {'path': empty_path, 'status': 'running'}, + } + deleted = list() + + with patch.object(crest_mod, 'check_job_status', lambda job_id: 'running'), \ + patch.object(crest_mod, 'delete_job', deleted.append): + crest_mod.monitor_crest_jobs(crest_jobs, check_interval=0, max_time_seconds=0) + + self.assertEqual(sorted(deleted), ['1', '2']) + self.assertEqual(crest_jobs['1']['status'], 'done') + self.assertEqual(crest_jobs['2']['status'], 'errored') + + references = {finished_path: {'xyz': OH_OH_XYZ, 'constraints': OH_OH_CONSTRAINTS}} + self.assertEqual(crest_mod.process_completed_jobs(crest_jobs, references), [OH_OH_XYZ]) + + def test_monitor_logs_a_cancellation_failure_and_keeps_going(self): + """One failed cancellation must not stop the remaining jobs from being cancelled.""" + crest_jobs = {'1': {'path': os.path.join(self.tmpdir.name, 'crest_0'), 'status': 'running'}, + '2': {'path': os.path.join(self.tmpdir.name, 'crest_1'), 'status': 'running'}, + } + deleted = list() + + def fake_delete(job_id): + deleted.append(job_id) + if job_id == '1': + raise RuntimeError(f'Could not delete job {job_id}') + + with patch.object(crest_mod, 'check_job_status', lambda job_id: 'running'), \ + patch.object(crest_mod, 'delete_job', fake_delete), \ + self.assertLogs('arc', level='ERROR') as log: + crest_mod.monitor_crest_jobs(crest_jobs, check_interval=0, max_time_seconds=0) + + self.assertEqual(sorted(deleted), ['1', '2']) + self.assertEqual(crest_jobs['1']['status'], 'errored') + self.assertEqual(crest_jobs['2']['status'], 'errored') + self.assertTrue(any('Could not cancel the timed-out CREST job 1' in message for message in log.output)) + + def test_failed_submission_does_not_collapse_jobs(self): + """submit_job() signals a failed submission as (None, None) or as ('errored', '').""" + submissions = [(None, None), ('errored', ''), ('running', '17'), ('running', '18')] + paths = [os.path.join(self.tmpdir.name, f'crest_{i}') for i in range(4)] + calls = {'n': 0} + + def fake_submit(path): + calls['n'] += 1 + return submissions[calls['n'] - 1] + + with patch.object(crest_mod, 'submit_job', fake_submit): + jobs = crest_mod.submit_crest_jobs(paths) + + self.assertNotIn(None, jobs) + self.assertNotIn('', jobs) + self.assertEqual(len(jobs), 2) + self.assertEqual(sorted(info['path'] for info in jobs.values()), sorted(paths[2:])) + + def test_process_completed_jobs_rejects_bad_geometries(self): + """Colliding, non-finite and empty geometries are rejected; a valid one is accepted.""" + crest_path = os.path.join(self.tmpdir.name, 'crest_0') + os.makedirs(crest_path) + crest_best_path = os.path.join(crest_path, 'crest_best.xyz') + jobs = {'123': {'path': crest_path, 'status': 'done'}} + references = {crest_path: {'xyz': OH_OH_XYZ, 'constraints': OH_OH_CONSTRAINTS}} + + def write(content): + with open(crest_best_path, 'w') as f: + f.write(content) + + colliding_coords = list(OH_OH_XYZ['coords']) + colliding_coords[3] = (colliding_coords[2][0] + 0.3, colliding_coords[2][1], colliding_coords[2][2]) + colliding_xyz = dict(OH_OH_XYZ) + colliding_xyz['coords'] = tuple(colliding_coords) + write(f"4\nCREST geometry\n{xyz_to_str(colliding_xyz)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), []) + + write("""4 +CREST geometry +O 0.00000000 -0.02752832 -1.20590500 +H 0.00000000 -0.02752832 -0.03383145 +O 0.00000000 -0.02752832 1.12142787 +H 0.00000000 nan 1.37454478 +""") + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), []) + + write('0\nCREST geometry\n') + with self.assertLogs('arc', level='WARNING') as log: + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), []) + self.assertTrue(any('is empty' in message for message in log.output), + msg=f'A zero-atom geometry must be rejected as empty:\n{log.output}') + + write(f"4\nCREST geometry\n{xyz_to_str(OH_OH_XYZ)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, references), [OH_OH_XYZ]) + + def test_process_completed_jobs_rejects_a_dissociated_reactive_triad(self): + """Do not accept a crest_best.xyz whose acceptor has separated from the transferring H.""" + dissociated_xyz = str_to_xyz("""O -1.1644 0.0000 0.0000 + H 0.0000 0.0000 0.0000 + O 4.9000 0.0000 0.0000 + H 5.8703 0.0000 0.0000""") + zeus_bad_xyz = str_to_xyz("""O -0.71236464 0.03765902 -0.02937463 + H -0.60136223 -0.77534746 0.43444583 + O 0.69301187 0.05895500 0.02917997 + H 0.90856791 -0.75830305 -0.43135588""") + crest_path = os.path.join(self.tmpdir.name, 'crest_0') + os.makedirs(crest_path) + crest_best_path = os.path.join(crest_path, 'crest_best.xyz') + with open(crest_best_path, 'w') as f: + f.write(f"4\nCREST geometry\n{xyz_to_str(dissociated_xyz)}\n") + + jobs = {'123': {'path': crest_path, 'status': 'done'}} + references = {crest_path: {'xyz': OH_OH_XYZ, 'constraints': OH_OH_CONSTRAINTS}} + self.assertEqual(crest_mod.process_completed_jobs(jobs, crest_references={}), []) + self.assertEqual(crest_mod.process_completed_jobs(jobs, crest_references=references), []) + + with open(crest_best_path, 'w') as f: + f.write(f"4\nCREST geometry\n{xyz_to_str(zeus_bad_xyz)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, crest_references=references), []) + + with open(crest_best_path, 'w') as f: + f.write(f"4\nCREST geometry\n{xyz_to_str(OH_OH_XYZ)}\n") + self.assertEqual(crest_mod.process_completed_jobs(jobs, crest_references=references), [OH_OH_XYZ]) + + def test_preserves_reactive_constraints_validates_indices(self): + """Malformed constraint indices are rejected rather than trusted.""" + valid = {'atoms': (0, 1, 2), 'distance_pairs': ((0, 1), (1, 2))} + self.assertTrue(crest_mod._preserves_reactive_constraints(xyz=OH_OH_XYZ, + reference_xyz=OH_OH_XYZ, + constraints=valid)) + for constraints in [{'atoms': (0, 1, 1), 'distance_pairs': ((0, 1),)}, + {'atoms': (0, 1, 9), 'distance_pairs': ((0, 1),)}, + {'atoms': (0, 1, -1), 'distance_pairs': ((0, 1),)}, + {'atoms': tuple(), 'distance_pairs': tuple()}, + {'atoms': (0, 1, 2), 'distance_pairs': ((0, 3),)}, + {'atoms': (0, 1, 2), 'distance_pairs': ((0, 1, 2),)}, + ]: + self.assertFalse(crest_mod._preserves_reactive_constraints(xyz=OH_OH_XYZ, + reference_xyz=OH_OH_XYZ, + constraints=constraints), + msg=f'{constraints} should have been rejected') + + +class TestCrestExecuteIncore(CrestTestCase): + """Tests for the execute_incore() orchestration.""" + + def setUp(self): + """Build a CH4 + OH reaction and a CREST adapter writing into the temporary directory.""" + super().setUp() + self.rxn = make_ch4_oh_reaction() + self.patch_crest_module() + self.adapter = CrestAdapter(job_type='tsg', + reactions=[self.rxn], + testing=True, + project='test_CrestAdapter', + project_directory=self.tmpdir.name, + ) + os.makedirs(self.adapter.local_path, exist_ok=True) + + def run_execute_incore(self, crest_results, seeds=CH4_OH_SEED, submitted=None): + """ + Run execute_incore() with the submission, monitoring and ingestion steps patched out. + + Args: + crest_results (list): The geometries process_completed_jobs() should return. + seeds: The seed entry, or list of seed entries, get_ts_seeds() should return. + submitted (list, optional): A list collecting the CREST job directories submitted. + + Returns: + list: The CREST job directories that were prepared. + """ + prepared = submitted if submitted is not None else list() + seed_list = seeds if isinstance(seeds, list) else [seeds] + + def fake_submit(crest_paths): + prepared.extend(crest_paths) + return {str(i): {'path': path, 'status': 'done'} for i, path in enumerate(crest_paths)} + + with patch.object(crest_mod, 'crest_available', lambda: True), \ + patch.object(crest_mod, 'get_ts_seeds', lambda **kwargs: list(seed_list)), \ + patch.object(crest_mod, 'submit_crest_jobs', fake_submit), \ + patch.object(crest_mod, 'monitor_crest_jobs', lambda *args, **kwargs: None), \ + patch.object(crest_mod, 'process_completed_jobs', + lambda *args, **kwargs: list(crest_results)): + self.adapter.execute_incore() + return prepared + + def test_execute_incore_appends_a_crest_ts_guess(self): + """A CREST geometry becomes a successful TSGuess and is saved to the job directory.""" + self.run_execute_incore(crest_results=[CH4_OH_RESULT_XYZ]) + + self.assertIsNotNone(self.rxn.ts_species) + self.assertEqual(len(self.rxn.ts_species.ts_guesses), 1) + ts_guess = self.rxn.ts_species.ts_guesses[0] + self.assertEqual(ts_guess.method.lower(), 'crest') + self.assertTrue(ts_guess.success) + self.assertEqual(ts_guess.family, 'H_Abstraction') + self.assertEqual(ts_guess.initial_xyz['symbols'], CH4_OH_RESULT_XYZ['symbols']) + for guess_coords, expected_coords in zip(ts_guess.initial_xyz['coords'], CH4_OH_RESULT_XYZ['coords']): + for guess_coord, expected_coord in zip(guess_coords, expected_coords): + self.assertAlmostEqual(guess_coord, expected_coord, places=6) + self.assertTrue(os.path.isfile(os.path.join(self.adapter.local_path, 'CREST_0.xyz'))) + + def test_execute_incore_merges_a_duplicate_geometry(self): + """A CREST geometry equal to an existing guess adds 'crest' instead of a new guess.""" + self.rxn.ts_species = ARCSpecies(label='TS', is_ts=True, charge=0, multiplicity=2) + existing = TSGuess(method='autotst', success=True, family='H_Abstraction', xyz=CH4_OH_RESULT_XYZ) + self.rxn.ts_species.append_ts_guess(existing) + + self.run_execute_incore(crest_results=[CH4_OH_RESULT_XYZ]) + + self.assertEqual(len(self.rxn.ts_species.ts_guesses), 1) + merged = self.rxn.ts_species.ts_guesses[0] + self.assertEqual(merged.method, 'autotst') + self.assertIn('crest', merged.method_sources) + + def test_execute_incore_caps_the_number_of_seeds(self): + """Only MAX_CREST_SEEDS seeds are submitted, and dropping the rest is logged.""" + seeds = list() + for index in range(crest_mod.MAX_CREST_SEEDS * 3): + coords = tuple(tuple(coord + 0.05 * index for coord in atom_coords) + for atom_coords in CH4_OH_SEED_XYZ['coords']) + xyz = {'symbols': CH4_OH_SEED_XYZ['symbols'], + 'isotopes': CH4_OH_SEED_XYZ['isotopes'], + 'coords': coords, + } + seeds.append(dict(CH4_OH_SEED, xyz=xyz)) + + submitted = list() + with self.assertLogs('arc', level='INFO') as log: + self.run_execute_incore(crest_results=[], seeds=seeds, submitted=submitted) + + self.assertEqual(len(submitted), crest_mod.MAX_CREST_SEEDS) + self.assertEqual(len(set(submitted)), crest_mod.MAX_CREST_SEEDS) + self.assertTrue(any(f'dropping the other {len(seeds) - crest_mod.MAX_CREST_SEEDS}' in message + for message in log.output), + msg=f'The dropped seeds were not logged:\n{log.output}') + + def test_execute_incore_skips_only_the_seed_that_could_not_be_prepared(self): + """A seed whose directory could not be prepared is dropped on its own. + + crest_ts_conformer_search() returns None for a directory holding a stale crest_best.xyz it + could not remove, which is a property of that one directory and says nothing about the + remaining seeds. + """ + seeds = list() + for index in range(2): + coords = tuple(tuple(coord + 0.5 * index for coord in atom_coords) + for atom_coords in CH4_OH_SEED_XYZ['coords']) + seeds.append(dict(CH4_OH_SEED, xyz={'symbols': CH4_OH_SEED_XYZ['symbols'], + 'isotopes': CH4_OH_SEED_XYZ['isotopes'], + 'coords': coords, + })) + calls = {'n': 0} + real_search = crest_mod.crest_ts_conformer_search + + def fake_search(*args, **kwargs): + calls['n'] += 1 + if calls['n'] == 1: + return None + return real_search(*args, **kwargs) + + submitted = list() + with patch.object(crest_mod, 'crest_ts_conformer_search', fake_search): + self.run_execute_incore(crest_results=[], seeds=seeds, submitted=submitted) + + self.assertEqual(calls['n'], 2) + self.assertEqual([os.path.basename(path) for path in submitted], ['crest_0_1']) + + def test_execute_incore_names_job_directories_per_reaction(self): + """The CREST job directory name carries the reaction index, so reactions cannot collide.""" + submitted = list() + self.run_execute_incore(crest_results=[], submitted=submitted) + self.assertEqual([os.path.basename(path) for path in submitted], ['crest_0_0']) + + def test_execute_incore_falls_back_to_an_external_ts_guess(self): + """With no heuristic seed, CREST is seeded from another adapter's successful TS guess.""" + self.rxn.ts_species = ARCSpecies(label='TS', is_ts=True, charge=0, multiplicity=2) + self.rxn.ts_species.append_ts_guess( + TSGuess(method='autotst', success=True, family='H_Abstraction', xyz=CH4_OH_SEED_XYZ) + ) + recorded = list() + + def fake_search(xyz_guess, **kwargs): + recorded.append({'xyz': xyz_guess, 'constraints': kwargs.get('constraints')}) + crest_path = os.path.join(self.tmpdir.name, f"crest_{kwargs.get('xyz_crest_int')}") + os.makedirs(crest_path, exist_ok=True) + return crest_path + + with patch.object(crest_mod, 'crest_available', lambda: True), \ + patch.object(crest_mod, 'get_ts_seeds', lambda **kwargs: []), \ + patch.object(crest_mod, 'crest_ts_conformer_search', fake_search), \ + patch.object(crest_mod, 'submit_crest_jobs', lambda paths: {}), \ + patch.object(crest_mod, 'monitor_crest_jobs', lambda *args, **kwargs: None), \ + patch.object(crest_mod, 'process_completed_jobs', lambda *args, **kwargs: []): + self.adapter.execute_incore() + + self.assertEqual(len(recorded), 1) + self.assertEqual(recorded[0]['xyz']['symbols'], CH4_OH_SEED_XYZ['symbols']) + constraints = recorded[0]['constraints'] + self.assertIsNotNone(constraints) + self.assertEqual(set(constraints['atoms']), {0, 1, 5}) + self.assertEqual(constraints['angle_atoms'][1], 1) + for pair in constraints['distance_pairs']: + self.assertEqual(len(pair), 2) + self.assertTrue(all(atom in constraints['atoms'] for atom in pair)) + + +class TestCrestHydrolysisConstraints(CrestTestCase): + """Tests for the four-center hydrolysis CREST constraint specification.""" + + def test_six_distance_pairs_for_every_hydrolysis_family(self): + """Each hydrolysis family pins all six pairwise distances among a, b, o and h1.""" + for family in HYDROLYSIS_FAMILIES: + rxn = types.SimpleNamespace(family=family) + constraints = get_wrapper_constraints(wrapper='crest', + reaction=rxn, + seed=make_hydrolysis_seed(family=family), + ) + self.assertIsNotNone(constraints, msg=f'No CREST constraints were built for {family}.') + self.assertEqual(constraints['atoms'], (2, 1, 8, 9)) + self.assertEqual(len(constraints['distance_pairs']), 6) + self.assertEqual({frozenset(pair) for pair in constraints['distance_pairs']}, + {frozenset(pair) for pair in HYDROLYSIS_DISTANCE_PAIRS}) + self.assertNotIn('angle_atoms', constraints) + + def test_six_pairs_are_derived_from_the_positional_metadata_alone(self): + """A seed carrying only the positional indices still yields the full six-distance core.""" + seed = make_hydrolysis_seed(reactive_atoms=False) + self.assertNotIn('reactive_atoms', seed['metadata']) + constraints = get_wrapper_constraints(wrapper='crest', + reaction=types.SimpleNamespace(family='carbonyl_based_hydrolysis'), + seed=seed, + ) + self.assertEqual(constraints['atoms'], (2, 1, 8, 9)) + self.assertEqual({frozenset(pair) for pair in constraints['distance_pairs']}, + {frozenset(pair) for pair in HYDROLYSIS_DISTANCE_PAIRS}) + + def test_seed_metadata_gains_the_named_reactive_atom_dict(self): + """get_ts_seeds() converts the positional hydrolysis indices into a named role dict.""" + rxn = types.SimpleNamespace(family='ether_hydrolysis') + with patch('arc.job.adapters.ts.heuristics.hydrolysis', + return_value=([HYDROLYSIS_SEED_XYZ], + ['ether_hydrolysis'], + [list(HYDROLYSIS_POSITIONAL_INDICES)])): + seeds = get_ts_seeds(reaction=rxn, base_adapter='heuristics') + self.assertEqual(len(seeds), 1) + self.assertEqual(seeds[0]['metadata']['reactive_atoms'], HYDROLYSIS_REACTIVE_ATOMS) + self.assertEqual(seeds[0]['metadata']['indices'], HYDROLYSIS_POSITIONAL_INDICES) + + def test_seed_metadata_accepts_an_already_named_index_mapping(self): + """A hydrolysis generator emitting a role-keyed mapping is honoured as-is.""" + named = {'a': 2, 'b': 1, 'e': 3, 'o': 8, 'h1': 9, 'd': 7} + rxn = types.SimpleNamespace(family='nitrile_hydrolysis') + with patch('arc.job.adapters.ts.heuristics.hydrolysis', + return_value=([HYDROLYSIS_SEED_XYZ], ['nitrile_hydrolysis'], [named])): + seeds = get_ts_seeds(reaction=rxn, base_adapter='heuristics') + self.assertEqual(seeds[0]['metadata']['reactive_atoms'], HYDROLYSIS_REACTIVE_ATOMS) + constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seeds[0]) + self.assertEqual(constraints['atoms'], (2, 1, 8, 9)) + + def test_constraints_are_none_for_missing_or_malformed_metadata(self): + """A hydrolysis seed whose roles cannot be resolved yields no constraints, with a warning.""" + rxn = types.SimpleNamespace(family='carbonyl_based_hydrolysis') + broken_metadata = [None, + {}, + {'indices': None}, + {'indices': [2, 1, 8]}, + {'indices': {'a': 2, 'b': 1, 'o': 8}}, + {'indices': [2, 1, 3, 8, 99, 7]}, + {'indices': [2, 1, 3, 8, 8, 7]}, + {'indices': [2, 1, 3, 9, 8, 7]}, + {'reactive_atoms': {'a': 2, 'b': 1, 'o': 9, 'h1': 8}}, + {'reactive_atoms': {'A': 2, 'H': 9, 'B': 1}}, + ] + for metadata in broken_metadata: + seed = {'xyz': HYDROLYSIS_SEED_XYZ, + 'family': 'carbonyl_based_hydrolysis', + 'metadata': metadata, + } + with self.assertLogs('arc', level='WARNING'): + constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) + self.assertIsNone(constraints, msg=f'Metadata {metadata!r} should not yield constraints.') + + def test_constraints_inp_holds_one_based_indices(self): + """A hydrolysis seed round trips into a $constrain block of six one-based distances.""" + self.patch_crest_module() + constraints = get_wrapper_constraints(wrapper='crest', + reaction=types.SimpleNamespace(family='carbonyl_based_hydrolysis'), + seed=make_hydrolysis_seed(), + ) + crest_path = crest_mod.crest_ts_conformer_search(xyz_guess=HYDROLYSIS_SEED_XYZ, + constraints=constraints, + path=self.tmpdir.name, + xyz_crest_int='hydro', + ) + with open(os.path.join(crest_path, 'constraints.inp')) as f: + constraints_text = f.read() + + self.assertIn('atoms: 3, 2, 9, 10\n', constraints_text) + for atom_1, atom_2 in HYDROLYSIS_DISTANCE_PAIRS: + self.assertIn(f'distance: {atom_1 + 1}, {atom_2 + 1}, auto', constraints_text, + msg=f'The {atom_1}--{atom_2} distance is missing from:\n{constraints_text}') + self.assertEqual(constraints_text.count('distance:'), 6) + self.assertNotIn('angle:', constraints_text) + self.assertIn('$metadyn\n atoms: 1, 4, 5, 6, 7, 8, 11\n', constraints_text) + self.assertNotIn('distance: 2, 1,', constraints_text) + self.assertNotIn('distance: 8, 9,', constraints_text) + + def test_a_hydrolysis_geometry_survives_the_reactive_constraint_check(self): + """The six-distance core accepts the seed geometry and rejects a dissociated one.""" + constraints = get_wrapper_constraints(wrapper='crest', + reaction=types.SimpleNamespace(family='nitrile_hydrolysis'), + seed=make_hydrolysis_seed(family='nitrile_hydrolysis'), + ) + self.assertTrue(crest_mod._preserves_reactive_constraints(xyz=HYDROLYSIS_SEED_XYZ, + reference_xyz=HYDROLYSIS_SEED_XYZ, + constraints=constraints, + )) + dissociated_coords = list(HYDROLYSIS_SEED_XYZ['coords']) + dissociated_coords[8] = (9.0, -0.06034836, 1.24530570) + dissociated_xyz = dict(HYDROLYSIS_SEED_XYZ, coords=tuple(dissociated_coords)) + self.assertFalse(crest_mod._preserves_reactive_constraints(xyz=dissociated_xyz, + reference_xyz=HYDROLYSIS_SEED_XYZ, + constraints=constraints, + )) + + +class TestCrestHydrolysisExecution(CrestTestCase): + """Tests for driving execute_incore() with hydrolysis seeds.""" + + def setUp(self): + """Build a methyl formate hydrolysis reaction and a CREST adapter.""" + super().setUp() + self.rxn = make_methyl_formate_hydrolysis_reaction() + self.patch_crest_module() + self.adapter = CrestAdapter(job_type='tsg', + reactions=[self.rxn], + testing=True, + project='test_CrestAdapter', + project_directory=self.tmpdir.name, + ) + os.makedirs(self.adapter.local_path, exist_ok=True) + + def run_with_seeds(self, seeds): + """ + Run execute_incore() over ``seeds`` and return the recorded constraint specifications. + + Args: + seeds (list): The seed entries get_ts_seeds() should return. + + Returns: + list: One entry per prepared CREST job, holding its constraint specification. + """ + recorded = list() + + def fake_search(xyz_guess, **kwargs): + recorded.append(kwargs.get('constraints')) + crest_path = os.path.join(self.tmpdir.name, f"crest_{kwargs.get('xyz_crest_int')}") + os.makedirs(crest_path, exist_ok=True) + return crest_path + + with patch.object(crest_mod, 'crest_available', lambda: True), \ + patch.object(crest_mod, 'get_ts_seeds', lambda **kwargs: list(seeds)), \ + patch.object(crest_mod, 'get_backup_ts_seeds', lambda *args, **kwargs: []), \ + patch.object(crest_mod, 'crest_ts_conformer_search', fake_search), \ + patch.object(crest_mod, 'submit_crest_jobs', lambda paths: {}), \ + patch.object(crest_mod, 'monitor_crest_jobs', lambda *args, **kwargs: None), \ + patch.object(crest_mod, 'process_completed_jobs', lambda *args, **kwargs: []): + self.adapter.execute_incore() + return recorded + + def test_a_hydrolysis_seed_reaches_crest_with_six_distances(self): + """CREST accepts the registered hydrolysis family and constrains the four-center core.""" + recorded = self.run_with_seeds([make_hydrolysis_seed()]) + self.assertEqual(len(recorded), 1) + self.assertEqual(recorded[0]['atoms'], (2, 1, 8, 9)) + self.assertEqual({frozenset(pair) for pair in recorded[0]['distance_pairs']}, + {frozenset(pair) for pair in HYDROLYSIS_DISTANCE_PAIRS}) + + def test_a_malformed_hydrolysis_seed_is_skipped_cleanly(self): + """A seed with unusable metadata is skipped and logged rather than aborting the run.""" + seed = make_hydrolysis_seed(indices=[2, 1, 8], reactive_atoms=False) + with self.assertLogs('arc', level='WARNING') as log: + recorded = self.run_with_seeds([seed]) + self.assertEqual(recorded, []) + self.assertTrue(any('CREST hydrolysis reactive atoms' in message for message in log.output), + msg=f'The unusable hydrolysis metadata was not logged:\n{log.output}') + self.assertTrue(any('Skipping this CREST seed' in message for message in log.output), + msg=f'The skipped seed was not logged:\n{log.output}') + + +class TestCrestHydrolysisRegistration(unittest.TestCase): + """Tests for the CREST family registration of the hydrolysis families.""" + + def test_crest_registered_for_every_hydrolysis_family(self): + """CREST is an eligible TS adapter for all three hydrolysis families.""" + for family in HYDROLYSIS_FAMILIES: + self.assertIn('crest', ts_adapters_by_rmg_family[family]) + + def test_heuristics_stays_registered_for_every_hydrolysis_family(self): + """Registering CREST must not displace the heuristics adapter.""" + for family in HYDROLYSIS_FAMILIES: + self.assertIn('heuristics', ts_adapters_by_rmg_family[family]) + + def test_hydrolysis_families_match_family_sets(self): + """The hub's local family tuple must stay in step with heuristics.FAMILY_SETS. + + seed_hub holds its own copy so that resolving the families does not import heuristics. + This test is what makes that duplication safe. + """ + from arc.job.adapters.ts.heuristics import FAMILY_SETS + from arc.job.adapters.ts.seed_hub import get_hydrolysis_families + + from_family_sets = sorted(family for families in FAMILY_SETS.values() for family in families) + self.assertEqual(sorted(get_hydrolysis_families()), from_family_sets) + self.assertEqual(sorted(HYDROLYSIS_FAMILIES), from_family_sets) + + def test_crest_is_not_a_default_ts_adapter(self): + """CREST stays opt-in: it is absent from the default ts_adapters list. + + An adapter runs only if it is both registered for the family and listed in ts_adapters, + so registering the hydrolysis families does not change ARC's default behaviour. + """ + self.assertNotIn('crest', [adapter.lower() for adapter in ts_adapters]) + + +if __name__ == '__main__': + unittest.main() diff --git a/arc/job/adapters/ts/heuristics.py b/arc/job/adapters/ts/heuristics.py index e0d227e56e..dcac20e4f4 100644 --- a/arc/job/adapters/ts/heuristics.py +++ b/arc/job/adapters/ts/heuristics.py @@ -21,31 +21,52 @@ import os from typing import TYPE_CHECKING, Any -from arc.common import (ARC_PATH, almost_equal_coords, get_angle_in_180_range, get_logger, is_angle_linear, - is_xyz_linear, key_by_val, read_yaml_file) +from arc.common import ( + ARC_PATH, + almost_equal_coords, + get_angle_in_180_range, + get_logger, + is_angle_linear, + is_xyz_linear, + key_by_val, + read_yaml_file, +) from arc.family import get_reaction_family_products from arc.job.adapter import JobAdapter from arc.job.adapters.common import _initialize_adapter, ts_adapters_by_rmg_family from arc.job.factory import register_job_adapter from arc.plotter import save_geo -from arc.species.converter import (compare_zmats, relocate_zmat_dummy_atoms_to_the_end, zmat_from_xyz, zmat_to_xyz, - add_atom_to_xyz_using_internal_coords, sorted_distances_of_atom) +from arc.species.converter import ( + add_atom_to_xyz_using_internal_coords, + compare_zmats, + relocate_zmat_dummy_atoms_to_the_end, + sorted_distances_of_atom, + zmat_from_xyz, + zmat_to_xyz, +) from arc.mapping.engine import map_two_species from arc.molecule.molecule import Molecule from arc.species.species import ARCSpecies, TSGuess, SpeciesError, colliding_atoms from arc.species.zmat import get_parameter_from_atom_indices, remove_zmat_atom_0, up_param, xyz_to_zmat from arc.species.vectors import calculate_angle +from arc.job.adapters.ts.seed_hub import get_ts_seeds if TYPE_CHECKING: from arc.level import Level from arc.reaction import ARCReaction - -FAMILY_SETS = {'hydrolysis_set_1': ['carbonyl_based_hydrolysis', 'ether_hydrolysis'], - 'hydrolysis_set_2': ['nitrile_hydrolysis']} +FAMILY_SETS = { + 'hydrolysis_set_1': ['carbonyl_based_hydrolysis', 'ether_hydrolysis'], + 'hydrolysis_set_2': ['nitrile_hydrolysis'], +} DIHEDRAL_INCREMENT = 30 +# Absolute reactive-core distances (Angstrom) for an H_Abstraction transition state A-H...B, +# taken as the medians over 1670 DFT-optimized H-abstraction transition states. +BREAKING_BOND_LENGTH = 1.29 +FORMING_BOND_LENGTH = 1.36 + ELECTRONEGATIVITIES = read_yaml_file(os.path.join(ARC_PATH, 'data', 'electronegativity.yml')) logger = get_logger() @@ -258,55 +279,60 @@ def execute_incore(self): multiplicity=rxn.multiplicity, ) - xyzs = list() - tsg, families = None, None - if rxn.family == 'H_Abstraction': - tsg = TSGuess(method='Heuristics') - tsg.tic() - xyzs = h_abstraction(reaction=rxn, dihedral_increment=self.dihedral_increment) - tsg.tok() - + tsg = TSGuess(method='Heuristics') + tsg.tic() + xyzs = get_ts_seeds( + reaction=rxn, + base_adapter='heuristics', + dihedral_increment=self.dihedral_increment, + ) + tsg.tok() if rxn.family in FAMILY_SETS['hydrolysis_set_1'] or rxn.family in FAMILY_SETS['hydrolysis_set_2']: - try: - tsg = TSGuess(method='Heuristics') - tsg.tic() - xyzs, families, indices = hydrolysis(reaction=rxn) - tsg.tok() - if not xyzs: - logger.warning(f'Heuristics TS search failed to generate any valid TS guesses for {rxn.label}.') - continue - except ValueError: + if not xyzs: + logger.warning( + f'Heuristics TS search failed to generate any valid TS guesses for {rxn.label}.' + ) continue - for method_index, xyz in enumerate(xyzs): + for method_index, xyz_entry in enumerate(xyzs): + xyz = xyz_entry.get("xyz") + method_label = xyz_entry.get("method", "Heuristics") + family = xyz_entry.get("family", rxn.family) + if xyz is None: + continue unique = True for other_tsg in rxn.ts_species.ts_guesses: if almost_equal_coords(xyz, other_tsg.initial_xyz): - if 'heuristics' not in other_tsg.method.lower(): - other_tsg.method += ' and Heuristics' + existing_sources = getattr(other_tsg, "method_sources", None) + if existing_sources is not None: + combined_sources = list(existing_sources) + [method_label] + else: + combined_sources = [other_tsg.method, method_label] + other_tsg.method_sources = TSGuess._normalize_method_sources(combined_sources) unique = False break if unique: - ts_guess = TSGuess(method='Heuristics', + ts_guess = TSGuess(method=method_label, method_index=method_index, t0=tsg.t0, execution_time=tsg.execution_time, success=True, - family=rxn.family if families is None else families[method_index], + family=family, xyz=xyz, ) rxn.ts_species.append_ts_guess(ts_guess) save_geo(xyz=xyz, path=self.local_path, - filename=f'Heuristics_{method_index}', + filename=f'{method_label}_{method_index}', format_='xyz', - comment=f'Heuristics {method_index}, family: {rxn.family}', + comment=f'{method_label} {method_index}, family: {rxn.family}', ) if len(self.reactions) < 5: - successes = len([tsg for tsg in rxn.ts_species.ts_guesses if tsg.success and 'heuristics' in tsg.method]) + successes = [tsg for tsg in rxn.ts_species.ts_guesses + if tsg.success and 'heuristics' in tsg.method.lower()] if successes: - logger.info(f'Heuristics successfully found {successes} TS guesses for {rxn.label}.') + logger.info(f'Heuristics successfully found {len(successes)} TS guesses for {rxn.label}.') else: logger.info(f'Heuristics did not find any successful TS guesses for {rxn.label}.') @@ -329,8 +355,8 @@ def combine_coordinates_with_redundant_atoms(xyz_1: dict[str, Any], h2: int, c: int | None = None, d: int | None = None, - r1_stretch: float = 1.2, - r2_stretch: float = 1.2, + r1_stretch: float | None = None, + r2_stretch: float | None = None, a2: float = 180.0, d2: float | None = None, d3: float | None = None, @@ -358,13 +384,13 @@ def combine_coordinates_with_redundant_atoms(xyz_1: dict[str, Any], ('R_1_0', None, None), # 1, atom C ('R_2_1', 'A_2_1_0', None), # 2 ('R_3_2', 'A_3_2_0', 'D_3_2_0_1'), # 3 - ('R_4_3' * r1_stretch, 'A_4_3_2', 'D_4_3_2_1')), # 4, H1 + ('R_4_3', 'A_4_3_2', 'D_4_3_2_1')), # 4, H1, R_4_3 is the breaking bond 'vars': {...}, 'map': {...}} zmat2 = {'symbols': ('H', 'H', 'H', 'H', 'C'), 'coords': ((None, None, None), # H2, redundant H atom, will be united with H1 - ('R_5_4' * r2_stretch, a2 (B-H-A) = 'A_5_4_0', d2 (B-H-A-C) = 'D_5_4_0_1'), # 5, atom B + ('R_5_4', a2 (B-H-A) = 'A_5_4_0', d2 (B-H-A-C) = 'D_5_4_0_1'), # 5, B, forming bond ('R_6_5', 'A_6_5_4', d3 (D-B-H-A) = 'D_6_5_4_0'), # 6, atom D ('R_7_6', 'A_7_6_4', 'D_7_6_4_5'), # 7 ('R_8_7', 'A_8_7_6', 'D_8_7_6_5')), # 8 @@ -384,10 +410,12 @@ def combine_coordinates_with_redundant_atoms(xyz_1: dict[str, Any], (atom C). d (int | None): The 0-index of an atom in ``xyz2`` connected to either B or H2 which is neither B nor H2 (atom D). - r1_stretch (float, optional): The factor by which to multiply (stretch/shrink) the bond length to the terminal - atom ``h1`` in ``xyz1`` (bond A-H1). - r2_stretch (float, optional): The factor by which to multiply (stretch/shrink) the bond length to the terminal - atom ``h2`` in ``xyz2`` (bond B-H2). + r1_stretch (float | None, optional): The factor by which to multiply (stretch/shrink) the bond length to the + terminal atom ``h1`` in ``xyz1`` (bond A-H1). If ``None``, the bond is set + to the absolute ``BREAKING_BOND_LENGTH`` instead of being multiplied. + r2_stretch (float | None, optional): The factor by which to multiply (stretch/shrink) the bond length to the + terminal atom ``h2`` in ``xyz2`` (bond B-H2). If ``None``, the bond is set + to the absolute ``FORMING_BOND_LENGTH`` instead of being multiplied. a2 (float, optional): The angle (in degrees) in the combined structure between atoms B-H-A (angle B-H-A). d2 (float | None): The dihedral angle (in degrees) between atoms B-H-A-C (dihedral B-H-A-C). This argument must be given only if the a2 angle is not linear, @@ -408,9 +436,15 @@ def combine_coordinates_with_redundant_atoms(xyz_1: dict[str, Any], is_a2_linear, is_mol_1_linear, a, b = _validate_combine_coordinates_with_redundant_atoms_args( xyz_1, xyz_2, mol_1, mol_2, h1, h2, a2, d2, d3, c, d) zmat_1, zmat_2 = generate_the_two_constrained_zmats(xyz_1, xyz_2, mol_1, mol_2, h1, h2, a, b, c, d) - # Stretch the A--H1 and B--H2 bonds. - stretch_zmat_bond(zmat=zmat_1, indices=(h1, a), stretch=r1_stretch) - stretch_zmat_bond(zmat=zmat_2, indices=(b, h2), stretch=r2_stretch) + # Set the breaking A--H1 and forming B--H2 bonds. + if r1_stretch is None: + set_zmat_bond(zmat=zmat_1, indices=(h1, a), length=BREAKING_BOND_LENGTH) + else: + stretch_zmat_bond(zmat=zmat_1, indices=(h1, a), stretch=r1_stretch) + if r2_stretch is None: + set_zmat_bond(zmat=zmat_2, indices=(b, h2), length=FORMING_BOND_LENGTH) + else: + stretch_zmat_bond(zmat=zmat_2, indices=(b, h2), stretch=r2_stretch) add_dummy = is_a2_linear and len(zmat_1['symbols']) > 2 and not is_mol_1_linear glue_params = determine_glue_params(zmat=zmat_1, add_dummy=add_dummy, @@ -552,6 +586,21 @@ def stretch_zmat_bond(zmat: dict, zmat['vars'][param] *= stretch +def set_zmat_bond(zmat: dict, + indices: tuple[int, int], + length: float): + """ + Set a bond in a zmat to an absolute length. + + Args: + zmat (dict): The zmat to process. + indices (tuple): A length 2 tuple with the 0-indices of the xyz (not zmat) atoms representing the bond to set. + length (float): The bond length to set, in Angstrom. + """ + param = get_parameter_from_atom_indices(zmat=zmat, indices=indices, xyz_indexed=True) + zmat['vars'][param] = length + + def determine_glue_params(zmat: dict, add_dummy: bool, h1: int, @@ -853,8 +902,8 @@ def are_h_abs_wells_reversed(rxn: ARCReaction, def h_abstraction(reaction: ARCReaction, - r1_stretch: float = 1.2, - r2_stretch: float = 1.2, + r1_stretch: float | None = None, + r2_stretch: float | None = None, a2: float = 180, dihedral_increment: int | None = None, ) -> list[dict]: @@ -863,10 +912,14 @@ def h_abstraction(reaction: ARCReaction, Args: reaction: An ARCReaction instance. - r1_stretch (float, optional): The factor by which to multiply (stretch/shrink) the bond length to the terminal - atom ``h1`` in ``xyz1`` (bond A-H1) relative to the respective well. - r2_stretch (float, optional): The factor by which to multiply (stretch/shrink) the bond length to the terminal - atom ``h2`` in ``xyz2`` (bond B-H2) relative to the respective well. + r1_stretch (float | None, optional): The factor by which to multiply (stretch/shrink) the bond length to the + terminal atom ``h1`` in ``xyz1`` (bond A-H1) relative to the respective + well. If ``None``, the bond is set to the absolute + ``BREAKING_BOND_LENGTH`` instead of being multiplied. + r2_stretch (float | None, optional): The factor by which to multiply (stretch/shrink) the bond length to the + terminal atom ``h2`` in ``xyz2`` (bond B-H2) relative to the respective + well. If ``None``, the bond is set to the absolute + ``FORMING_BOND_LENGTH`` instead of being multiplied. a2 (float, optional): The angle (in degrees) in the combined structure between atoms B-H-A (angle B-H-A). dihedral_increment (int, optional): The dihedral increment to use for B-H-A-C and D-B-H-C dihedral scans. @@ -882,6 +935,11 @@ def h_abstraction(reaction: ARCReaction, return xyz_guesses reactants_reversed, products_reversed = are_h_abs_wells_reversed(rxn=reaction, product_dict=reaction.product_dicts[0]) for product_dict in reaction.product_dicts: + reactive_atoms = { + 'A': product_dict['r_label_map']['*1'], + 'H': product_dict['r_label_map']['*2'], + 'B': product_dict['r_label_map']['*3'], + } # Identify R1H and R2H in the "R1H + R2 <=> R1 + R2H" or "R2 + R1H <=> R2H + R1" reaction # The expected RMG atom labels are: R(*1)-H(*2) + R(*3)j <=> R(*1)j + R(*3)-H(*2). # They appear in each product_dict under the 'r_label_map' key. @@ -955,7 +1013,12 @@ def h_abstraction(reaction: ARCReaction, else: # This TS is unique, and has no atom collisions. zmats.append(zmat_guess) - xyz_guesses.append(xyz_guess) + xyz_guesses.append({ + "xyz": xyz_guess, + "method": "Heuristics", + "metadata": {"reactive_atoms": reactive_atoms}, + }) + return xyz_guesses @@ -990,9 +1053,11 @@ def hydrolysis(reaction: ARCReaction) -> tuple[list[dict], list[dict], list[int] is_set_1 = reaction_family in hydrolysis_parameters["family_sets"]["set_1"] is_set_2 = reaction_family in hydrolysis_parameters["family_sets"]["set_2"] - main_reactant, water, initial_xyz, xyz_indices = extract_reactant_and_indices(reaction, - product_dict, - is_set_1) + main_reactant, water, initial_xyz, xyz_indices = extract_reactant_and_indices( + reaction, + product_dict, + is_set_1, + ) base_xyz_indices = { "a": xyz_indices["a"], "b": xyz_indices["b"], @@ -1002,9 +1067,19 @@ def hydrolysis(reaction: ARCReaction) -> tuple[list[dict], list[dict], list[int] } adjustments_to_try = [False, True] if dihedrals_to_change_num == 1 else [True] for adjust_dihedral in adjustments_to_try: - chosen_xyz_indices, xyz_guesses, zmats_total, n_dihedrals_found = process_chosen_d_indices(initial_xyz, base_xyz_indices, xyz_indices, - hydrolysis_parameters,reaction_family, water, zmats_total, is_set_1, is_set_2, - dihedrals_to_change_num, should_adjust_dihedral=adjust_dihedral) + chosen_xyz_indices, xyz_guesses, zmats_total, n_dihedrals_found = process_chosen_d_indices( + initial_xyz, + base_xyz_indices, + xyz_indices, + hydrolysis_parameters, + reaction_family, + water, + zmats_total, + is_set_1, + is_set_2, + dihedrals_to_change_num, + should_adjust_dihedral=adjust_dihedral, + ) max_dihedrals_found = max(max_dihedrals_found, n_dihedrals_found) if xyz_guesses: xyz_guesses_total.extend(xyz_guesses) @@ -1018,8 +1093,8 @@ def hydrolysis(reaction: ARCReaction) -> tuple[list[dict], list[dict], list[int] condition_met = len(xyz_guesses_total) > 0 nitrile_in_inputs = any( - (pd.get("family") == "nitrile_hydrolysis") or - (isinstance(pd.get("family"), list) and "nitrile_hydrolysis" in pd.get("family")) + (pd.get("family") == "nitrile_hydrolysis") + or (isinstance(pd.get("family"), list) and "nitrile_hydrolysis" in pd.get("family")) for pd in product_dicts ) nitrile_already_found = any(fam == "nitrile_hydrolysis" for fam in reaction_families) @@ -1035,9 +1110,11 @@ def hydrolysis(reaction: ARCReaction) -> tuple[list[dict], list[dict], list[int] is_set_1 = reaction_family in hydrolysis_parameters["family_sets"]["set_1"] is_set_2 = reaction_family in hydrolysis_parameters["family_sets"]["set_2"] - main_reactant, water, initial_xyz, xyz_indices = extract_reactant_and_indices(reaction, - product_dict, - is_set_1) + main_reactant, water, initial_xyz, xyz_indices = extract_reactant_and_indices( + reaction, + product_dict, + is_set_1, + ) base_xyz_indices = { "a": xyz_indices["a"], "b": xyz_indices["b"], @@ -1051,10 +1128,18 @@ def hydrolysis(reaction: ARCReaction) -> tuple[list[dict], list[dict], list[int] break dihedrals_to_change_num += 1 chosen_xyz_indices, xyz_guesses, zmats_total, n_dihedrals_found = process_chosen_d_indices( - initial_xyz, base_xyz_indices, xyz_indices, - hydrolysis_parameters, reaction_family, water, zmats_total, is_set_1, is_set_2, - dihedrals_to_change_num, should_adjust_dihedral=True, - allow_nitrile_dihedrals=True + initial_xyz, + base_xyz_indices, + xyz_indices, + hydrolysis_parameters, + reaction_family, + water, + zmats_total, + is_set_1, + is_set_2, + dihedrals_to_change_num, + should_adjust_dihedral=True, + allow_nitrile_dihedrals=True, ) max_dihedrals_found = max(max_dihedrals_found, n_dihedrals_found) @@ -1086,11 +1171,13 @@ def get_products_and_check_families(reaction: ARCReaction) -> tuple[list[dict], consider_arc_families=True, ) carbonyl_based_present = any( - "carbonyl_based_hydrolysis" in (d.get("family", []) if isinstance(d.get("family"), list) else [d.get("family")]) + "carbonyl_based_hydrolysis" + in (d.get("family", []) if isinstance(d.get("family"), list) else [d.get("family")]) for d in product_dicts ) ether_present = any( - "ether_hydrolysis" in (d.get("family", []) if isinstance(d.get("family"), list) else [d.get("family")]) + "ether_hydrolysis" + in (d.get("family", []) if isinstance(d.get("family"), list) else [d.get("family")]) for d in product_dicts ) @@ -1166,11 +1253,13 @@ def extract_reactant_and_indices(reaction: ARCReaction, main_reactant, a_xyz_index, b_xyz_index, - two_neighbors + two_neighbors, ) except ValueError as e: - raise ValueError(f"Failed to determine neighbors by electronegativity for atom {a_xyz_index} " - f"in species {main_reactant.label}: {e}") + raise ValueError( + f"Failed to determine neighbors by electronegativity for atom {a_xyz_index} " + f"in species {main_reactant.label}: {e}" + ) o_index = len(main_reactant.mol.atoms) h1_index = o_index + 1 @@ -1181,7 +1270,7 @@ def extract_reactant_and_indices(reaction: ARCReaction, "e": e_xyz_index, "d": d_xyz_indices, "o": o_index, - "h1": h1_index + "h1": h1_index, } return main_reactant, water, initial_xyz, xyz_indices @@ -1226,11 +1315,18 @@ def process_chosen_d_indices(initial_xyz: dict, """ max_dihedrals_found = 0 for d_index in xyz_indices.get("d", []) or [None]: - chosen_xyz_indices = {**base_xyz_indices, "d": d_index} if d_index is not None else {**base_xyz_indices, - "d": None} + chosen_xyz_indices = {**base_xyz_indices, "d": d_index} if d_index is not None else { + **base_xyz_indices, + "d": None, + } current_zmat, zmat_indices = setup_zmat_indices(initial_xyz, chosen_xyz_indices) - matches = get_matching_dihedrals(current_zmat, zmat_indices['a'], zmat_indices['b'], - zmat_indices['e'], zmat_indices['d']) + matches = get_matching_dihedrals( + current_zmat, + zmat_indices['a'], + zmat_indices['b'], + zmat_indices['e'], + zmat_indices['d'], + ) max_dihedrals_found = max(max_dihedrals_found, len(matches)) if should_adjust_dihedral and dihedrals_to_change_num > len(matches): continue @@ -1248,22 +1344,28 @@ def process_chosen_d_indices(initial_xyz: dict, zmat_variants = generate_dihedral_variants(current_zmat, indices, adjustment_factors) if zmat_variants: adjusted_zmats.extend(zmat_variants) - if not adjusted_zmats: - pass - else: + if adjusted_zmats: zmats_to_process = adjusted_zmats ts_guesses_list = [] for zmat_to_process in zmats_to_process: ts_guesses, updated_zmats = process_family_specific_adjustments( - is_set_1, is_set_2, reaction_family, hydrolysis_parameters, - zmat_to_process, water, chosen_xyz_indices, zmats_total) + is_set_1, + is_set_2, + reaction_family, + hydrolysis_parameters, + zmat_to_process, + water, + chosen_xyz_indices, + zmats_total, + ) zmats_total = updated_zmats ts_guesses_list.extend(ts_guesses) if attempted_dihedral_adjustments and not ts_guesses_list and ( - reaction_family != 'nitrile_hydrolysis' or allow_nitrile_dihedrals): - flipped_zmats= [] + reaction_family != 'nitrile_hydrolysis' or allow_nitrile_dihedrals + ): + flipped_zmats = [] adjustment_factors = [15, 25, 35, 45, 55] for indices in indices_list: flipped_variants = generate_dihedral_variants(current_zmat, indices, adjustment_factors, flip=True) @@ -1271,8 +1373,14 @@ def process_chosen_d_indices(initial_xyz: dict, for zmat_to_process in flipped_zmats: ts_guesses, updated_zmats = process_family_specific_adjustments( - is_set_1, is_set_2, reaction_family, hydrolysis_parameters, - zmat_to_process, water, chosen_xyz_indices, zmats_total + is_set_1, + is_set_2, + reaction_family, + hydrolysis_parameters, + zmat_to_process, + water, + chosen_xyz_indices, + zmats_total, ) zmats_total = updated_zmats ts_guesses_list.extend(ts_guesses) @@ -1342,8 +1450,11 @@ def get_neighbors_by_electronegativity(spc: ARCSpecies, Raises: ValueError: If the atom has no valid neighbors. """ - neighbors = [neighbor for neighbor in spc.mol.atoms[atom_index].edges.keys() - if spc.mol.atoms.index(neighbor) != exclude_index] + neighbors = [ + neighbor + for neighbor in spc.mol.atoms[atom_index].edges.keys() + if spc.mol.atoms.index(neighbor) != exclude_index + ] if not neighbors: raise ValueError(f"Atom at index {atom_index} has no valid neighbors.") @@ -1357,12 +1468,17 @@ def get_neighbor_total_electronegativity(neighbor: 'Atom') -> float: float: The total electronegativity of the neighbor """ return sum( - ELECTRONEGATIVITIES[n.symbol] * neighbor.edges[n].order - for n in neighbor.edges.keys() + ELECTRONEGATIVITIES[n.symbol] * neighbor.edges[n].order for n in neighbor.edges.keys() ) - effective_electronegativities = [(ELECTRONEGATIVITIES[n.symbol] * spc.mol.atoms[atom_index].edges[n].order, - get_neighbor_total_electronegativity(n), n ) for n in neighbors] + effective_electronegativities = [ + ( + ELECTRONEGATIVITIES[n.symbol] * spc.mol.atoms[atom_index].edges[n].order, + get_neighbor_total_electronegativity(n), + n, + ) + for n in neighbors + ] effective_electronegativities.sort(reverse=True, key=lambda x: (x[0], x[1])) sorted_neighbors = [spc.mol.atoms.index(n[2]) for n in effective_electronegativities] most_electronegative = sorted_neighbors[0] @@ -1389,7 +1505,7 @@ def setup_zmat_indices(initial_xyz: dict, 'a': key_by_val(initial_zmat.get('map', {}), xyz_indices['a']), 'b': key_by_val(initial_zmat.get('map', {}), xyz_indices['b']), 'e': key_by_val(initial_zmat.get('map', {}), xyz_indices['e']), - 'd': key_by_val(initial_zmat.get('map', {}), xyz_indices['d']) if xyz_indices['d'] is not None else None + 'd': key_by_val(initial_zmat.get('map', {}), xyz_indices['d']) if xyz_indices['d'] is not None else None, } return initial_zmat, zmat_indices @@ -1400,15 +1516,15 @@ def generate_dihedral_variants(zmat: dict, flip: bool = False, tolerance_degrees: float = 10.0) -> list[dict]: """ - Create variants of a Z-matrix by adjusting dihedral angles using multiple adjustment factors. + Create variants of a Z-matrix by adjusting dihedral angles using multiple adjustment factors. This function creates variants of the Z-matrix using different adjustment factors: - 1. Retrieve the current dihedral value and normalize it to the (-180°, 180°] range. - 2. For each adjustment factor, slightly push the angle away from 0° or ±180° to avoid - unstable, boundary configurations. - 3. If `flip=True`, the same procedure is applied starting from a flipped - (180°-shifted) baseline angle. - 4. Each adjusted or flipped variant is deep-copied to ensure independence. + 1. Retrieve the current dihedral value and normalize it to the (-180°, 180°] range. + 2. For each adjustment factor, slightly push the angle away from 0° or ±180° to avoid + unstable, boundary configurations. + 3. If `flip=True`, the same procedure is applied starting from a flipped + (180°-shifted) baseline angle. + 4. Each adjusted or flipped variant is deep-copied to ensure independence. Args: zmat (dict): The initial Z-matrix. @@ -1416,7 +1532,8 @@ def generate_dihedral_variants(zmat: dict, adjustment_factors (list[float], optional): List of factors to try. flip (bool, optional): Whether to start from a flipped (180°) baseline dihedral angle. Defaults to False. - tolerance_degrees (float, optional): Tolerance (in degrees) for detecting angles near 0° or ±180°. Defaults to 10.0. + tolerance_degrees (float, optional): Tolerance (in degrees) for detecting angles near 0° or ±180°. + Defaults to 10.0. Returns: list[dict]: List of Z-matrix variants with adjusted dihedral angles. @@ -1442,8 +1559,9 @@ def push_up_dihedral(val: float, adj_factor: float) -> float: seed_value = normalized_value if flip: seed_value = get_angle_in_180_range(normalized_value + 180.0) - boundary_like = ((abs(seed_value) < tolerance_degrees) - or (180 - tolerance_degrees <= abs(seed_value) <= 180+tolerance_degrees)) + boundary_like = (abs(seed_value) < tolerance_degrees) or ( + 180 - tolerance_degrees <= abs(seed_value) <= 180 + tolerance_degrees + ) if boundary_like: for factor in adjustment_factors: variant = copy.deepcopy(zmat) @@ -1486,11 +1604,13 @@ def get_matching_dihedrals(zmat: dict, return matches -def stretch_ab_bond(initial_zmat: 'dict', - xyz_indices: 'dict', - zmat_indices: 'dict', - hydrolysis_parameters: 'dict', - reaction_family: str) -> None: +def stretch_ab_bond( + initial_zmat: dict, + xyz_indices: dict, + zmat_indices: dict, + hydrolysis_parameters: dict, + reaction_family: str, +) -> None: """ Stretch the bond between atoms a and b in the Z-matrix based on the reaction family parameters. @@ -1530,7 +1650,7 @@ def process_family_specific_adjustments(is_set_1: bool, xyz_indices: dict, zmats_total: list[dict]) -> tuple[list[dict], list[dict]]: """ - Process specific adjustments for different hydrolysis reaction families if needed, then generate TS guesses . + Process specific adjustments for different hydrolysis reaction families if needed, then generate TS guesses. Args: is_set_1 (bool): Whether the reaction belongs to parameter set 1. @@ -1548,21 +1668,34 @@ def process_family_specific_adjustments(is_set_1: bool, Raises: ValueError: If the reaction family is not supported. """ - a_xyz, b_xyz, e_xyz, o_xyz, h1_xyz, d_xyz= xyz_indices.values() + a_xyz, b_xyz, e_xyz, o_xyz, h1_xyz, d_xyz = xyz_indices.values() r_atoms = [a_xyz, o_xyz, o_xyz] a_atoms = [[b_xyz, a_xyz], [a_xyz, o_xyz], [h1_xyz, o_xyz]] - d_atoms = ([[e_xyz, d_xyz, a_xyz], [b_xyz, a_xyz, o_xyz], [a_xyz, h1_xyz, o_xyz]] - if d_xyz is not None else - [[e_xyz, b_xyz, a_xyz], [b_xyz, a_xyz, o_xyz], [a_xyz, h1_xyz, o_xyz]]) + d_atoms = ( + [[e_xyz, d_xyz, a_xyz], [b_xyz, a_xyz, o_xyz], [a_xyz, h1_xyz, o_xyz]] + if d_xyz is not None + else [[e_xyz, b_xyz, a_xyz], [b_xyz, a_xyz, o_xyz], [a_xyz, h1_xyz, o_xyz]] + ) r_value = hydrolysis_parameters['family_parameters'][str(reaction_family)]['r_value'] a_value = hydrolysis_parameters['family_parameters'][str(reaction_family)]['a_value'] d_values = hydrolysis_parameters['family_parameters'][str(reaction_family)]['d_values'] if is_set_1 or is_set_2: initial_xyz = zmat_to_xyz(initial_zmat) - return generate_hydrolysis_ts_guess(initial_xyz, xyz_indices.values(), water, r_atoms, a_atoms, d_atoms, - r_value, a_value, d_values, zmats_total, is_set_1, - threshold=0.6 if reaction_family == 'nitrile_hydrolysis' else 0.8) + return generate_hydrolysis_ts_guess( + initial_xyz, + xyz_indices.values(), + water, + r_atoms, + a_atoms, + d_atoms, + r_value, + a_value, + d_values, + zmats_total, + is_set_1, + threshold=0.6 if reaction_family == 'nitrile_hydrolysis' else 0.8, + ) else: raise ValueError(f"Family {reaction_family} not supported for hydrolysis TS guess generation.") @@ -1602,7 +1735,7 @@ def generate_hydrolysis_ts_guess(initial_xyz: dict, """ xyz_guesses = [] - for index, d_value in enumerate(d_values): + for d_value in d_values: xyz_guess = copy.deepcopy(initial_xyz) for i in range(3): xyz_guess = add_atom_to_xyz_using_internal_coords( @@ -1613,18 +1746,18 @@ def generate_hydrolysis_ts_guess(initial_xyz: dict, d_indices=d_atoms[i], r_value=r_value[i], a_value=a_value[i], - d_value=d_value[i] + d_value=d_value[i], ) - a_xyz, b_xyz, e_xyz, o_xyz, h1_xyz, d_xyz= xyz_indices - are_valid_bonds=check_ts_bonds(xyz_guess, [o_xyz, h1_xyz, h1_xyz+1, a_xyz, b_xyz]) - colliding=colliding_atoms(xyz_guess, threshold=threshold) + a_xyz, b_xyz, e_xyz, o_xyz, h1_xyz, d_xyz = xyz_indices + are_valid_bonds = check_ts_bonds(xyz_guess, [o_xyz, h1_xyz, h1_xyz + 1, a_xyz, b_xyz]) + colliding = colliding_atoms(xyz_guess, threshold=threshold) duplicate = any(compare_zmats(existing, xyz_to_zmat(xyz_guess)) for existing in zmats_total) if is_set_1: - dihedral_edao=[e_xyz, d_xyz, a_xyz, o_xyz] - dao_is_linear=check_dao_angle(dihedral_edao, xyz_guess) + dihedral_edao = [e_xyz, d_xyz, a_xyz, o_xyz] + dao_is_linear = check_dao_angle(dihedral_edao, xyz_guess) else: - dao_is_linear=False + dao_is_linear = False if xyz_guess is not None and not colliding and not duplicate and are_valid_bonds and not dao_is_linear: xyz_guesses.append(xyz_guess) zmats_total.append(xyz_to_zmat(xyz_guess)) @@ -1645,7 +1778,7 @@ def check_dao_angle(d_indices: list[int], xyz_guess: dict) -> bool: """ angle_indices = [d_indices[1], d_indices[2], d_indices[3]] angle_value = calculate_angle(xyz_guess, angle_indices) - norm_value=(angle_value + 180) % 180 + norm_value = (angle_value + 180) % 180 return (norm_value < 10) or (norm_value > 170) @@ -1660,7 +1793,7 @@ def check_ts_bonds(transition_state_xyz: dict, tested_atom_indices: list) -> boo Returns: bool: Whether the transition state guess has the expected water-related bonds. """ - oxygen_index, h1_index, h2_index, a_index, b_index= tested_atom_indices + oxygen_index, h1_index, h2_index, a_index, b_index = tested_atom_indices oxygen_bonds = sorted_distances_of_atom(transition_state_xyz, oxygen_index) h1_bonds = sorted_distances_of_atom(transition_state_xyz, h1_index) h2_bonds = sorted_distances_of_atom(transition_state_xyz, h2_index) @@ -1679,10 +1812,12 @@ def check_oxygen_bonds(bonds): return rel_error <= 0.1 return False - oxygen_has_valid_bonds = (oxygen_bonds[0][0] == h2_index and check_oxygen_bonds(oxygen_bonds)) - h1_has_valid_bonds = (h1_bonds[0][0] in {oxygen_index, b_index}and h1_bonds[1][0] in {oxygen_index, b_index}) + oxygen_has_valid_bonds = oxygen_bonds[0][0] == h2_index and check_oxygen_bonds(oxygen_bonds) + h1_has_valid_bonds = (h1_bonds[0][0] in {oxygen_index, b_index}) and ( + h1_bonds[1][0] in {oxygen_index, b_index} + ) h2_has_valid_bonds = h2_bonds[0][0] == oxygen_index return oxygen_has_valid_bonds and h1_has_valid_bonds and h2_has_valid_bonds -register_job_adapter('heuristics', HeuristicsAdapter) +register_job_adapter("heuristics", HeuristicsAdapter) diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py index ea7de2d8b5..e138010f23 100644 --- a/arc/job/adapters/ts/heuristics_test.py +++ b/arc/job/adapters/ts/heuristics_test.py @@ -10,10 +10,14 @@ import os import shutil import unittest +from types import SimpleNamespace +from unittest.mock import patch from arc.common import ARC_TESTING_PATH, almost_equal_coords from arc.family import get_reaction_family_products -from arc.job.adapters.ts.heuristics import (HeuristicsAdapter, +from arc.job.adapters.ts.heuristics import (BREAKING_BOND_LENGTH, + FORMING_BOND_LENGTH, + HeuristicsAdapter, are_h_abs_wells_reversed, combine_coordinates_with_redundant_atoms, determine_glue_params, @@ -22,25 +26,41 @@ get_modified_params_from_zmat_2, get_new_map_based_on_zmat_1, get_new_zmat_2_map, + set_zmat_bond, stretch_zmat_bond, get_main_reactant_and_water_from_hydrolysis_reaction, setup_zmat_indices, get_neighbors_by_electronegativity, get_matching_dihedrals, generate_dihedral_variants, + h_abstraction, check_dao_angle, check_ts_bonds, h_abstraction, ) +from arc.job.adapters.ts.seed_hub import (H_ATOM_BOND_CUTOFF, + _get_h_abs_atoms_from_xyz, + _h_atom_bond_cutoff, + _is_valid_h_abs_atom_assignment, + get_ts_seeds, + get_wrapper_constraints, + ) from arc.reaction import ARCReaction from arc.species.converter import str_to_xyz, zmat_to_xyz, zmat_from_xyz -from arc.species.species import ARCSpecies +from arc.species.species import ARCSpecies, TSGuess +from arc.species.vectors import calculate_angle, calculate_distance from arc.species.zmat import _compare_zmats, get_parameter_from_atom_indices from arc.species.species import check_isomorphism from arc.species.zmat import remove_zmat_atom_0 from arc.species.converter import relocate_zmat_dummy_atoms_to_the_end +# The empirical H abstraction reactive-core distances (in Angstrom) that the heuristics adapter is expected to +# reproduce. These literals are intentionally not imported from the heuristics module, so that a change to the +# module constants is caught here. +EXPECTED_BREAKING_BOND_LENGTH = 1.29 +EXPECTED_FORMING_BOND_LENGTH = 1.36 + class TestHeuristicsAdapter(unittest.TestCase): """ @@ -457,6 +477,91 @@ def setUpClass(cls): 'map': {0: 7, 1: 2, 2: 1, 3: 3, 4: 'X9', 5: 0, 6: 'X10', 7: 4, 8: 'X11', 9: 5, 10: 6, 11: 'X12', 12: 8}} + def assert_h_abs_reactive_core(self, xyz, a, h, b, a2=180.0, msg=''): + """ + Assert that an H abstraction TS guess has the expected reactive core A-H...B. + + Args: + xyz (dict): The Cartesian coordinates of the TS guess. + a (int): The 0-index of the donor heavy atom A. + h (int): The 0-index of the migrating hydrogen atom H. + b (int): The 0-index of the acceptor heavy atom B. + a2 (float, optional): The expected A-H-B angle in degrees. + msg (str, optional): An identifying message reported on an assertion failure. + """ + self.assertEqual(xyz['symbols'][h], 'H', msg=msg) + self.assertAlmostEqual(calculate_distance(coords=xyz, atoms=[a, h]), + EXPECTED_BREAKING_BOND_LENGTH, places=3, msg=f'{msg} (breaking A-H bond)') + self.assertAlmostEqual(calculate_distance(coords=xyz, atoms=[h, b]), + EXPECTED_FORMING_BOND_LENGTH, places=3, msg=f'{msg} (forming H-B bond)') + self.assertAlmostEqual(calculate_angle(coords=xyz, atoms=[a, h, b]), a2, places=2, msg=f'{msg} (A-H-B angle)') + + def test_bond_length_constants(self): + """Test the empirical H abstraction reactive-core bond length constants.""" + self.assertAlmostEqual(BREAKING_BOND_LENGTH, EXPECTED_BREAKING_BOND_LENGTH, places=6) + self.assertAlmostEqual(FORMING_BOND_LENGTH, EXPECTED_FORMING_BOND_LENGTH, places=6) + self.assertLess(BREAKING_BOND_LENGTH, FORMING_BOND_LENGTH) + + def test_set_zmat_bond(self): + """Test the set_zmat_bond() function.""" + zmat2_copy = copy.deepcopy(self.zmat_2) + self.assertNotAlmostEqual(self.zmat_2['vars']['R_2_1'], 1.29, places=3) + set_zmat_bond(zmat=zmat2_copy, indices=(1, 0), length=1.29) + self.assertAlmostEqual(zmat2_copy['vars']['R_2_1'], 1.29, places=10) + set_zmat_bond(zmat=zmat2_copy, indices=(1, 0), length=1.36) + self.assertAlmostEqual(zmat2_copy['vars']['R_2_1'], 1.36, places=10) + + def test_h_abstraction_reactive_core_distances(self): + """Test that h_abstraction() places the reactive core at the empirical absolute distances.""" + ch4 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) + ch3 = ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz) + oh = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) + h2o = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) + c2h6 = ARCSpecies(label='C2H6', smiles='CC', xyz=self.c2h6_xyz) + c2h5 = ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz) + for label, rxn in [('CH4 + OH', ARCReaction(r_species=[ch4, oh], p_species=[ch3, h2o])), + ('C2H6 + OH', ARCReaction(r_species=[c2h6, oh], p_species=[c2h5, h2o]))]: + seeds = h_abstraction(reaction=rxn, dihedral_increment=120) + self.assertTrue(len(seeds)) + for i, seed in enumerate(seeds): + reactive_atoms = seed['metadata']['reactive_atoms'] + self.assert_h_abs_reactive_core(xyz=seed['xyz'], + a=reactive_atoms['A'], + h=reactive_atoms['H'], + b=reactive_atoms['B'], + msg=f'{label} seed {i}', + ) + + def test_h_abstraction_forming_bond_is_element_independent(self): + """Test that the forming H-B bond does not depend on the acceptor element.""" + ch4 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) + ch3 = ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz) + oh = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) + h2o = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) + rxn_o_acceptor = ARCReaction(r_species=[ch4, oh], p_species=[ch3, h2o]) + rxn_h_acceptor = ARCReaction(r_species=[ch4, self.h], p_species=[ch3, self.h2]) + forming_bonds = dict() + for label, rxn in [('O acceptor', rxn_o_acceptor), ('H acceptor', rxn_h_acceptor)]: + seed = h_abstraction(reaction=rxn, dihedral_increment=120)[0] + reactive_atoms = seed['metadata']['reactive_atoms'] + forming_bonds[label] = calculate_distance(coords=seed['xyz'], + atoms=[reactive_atoms['H'], reactive_atoms['B']]) + self.assert_h_abs_reactive_core(xyz=seed['xyz'], + a=reactive_atoms['A'], + h=reactive_atoms['H'], + b=reactive_atoms['B'], + msg=label, + ) + self.assertAlmostEqual(forming_bonds['O acceptor'], forming_bonds['H acceptor'], places=6) + + # An explicit multiplicative stretch remains available, and is element-dependent. + stretched_seed = h_abstraction(reaction=ARCReaction(r_species=[ch4, oh], p_species=[ch3, h2o]), + r1_stretch=1.2, r2_stretch=1.2, dihedral_increment=120)[0] + reactive_atoms = stretched_seed['metadata']['reactive_atoms'] + self.assertAlmostEqual(calculate_distance(coords=stretched_seed['xyz'], + atoms=[reactive_atoms['H'], reactive_atoms['B']]), + 1.1628, places=3) + def test_heuristics_for_h_abstraction_1(self): """ Test that ARC can generate TS guesses based on heuristics for H Abstraction reactions. @@ -476,10 +581,13 @@ def test_heuristics_for_h_abstraction_1(self): self.assertEqual(rxn1.ts_species.charge, 0) self.assertEqual(rxn1.ts_species.multiplicity, 3) self.assertEqual(len(rxn1.ts_species.ts_guesses), 2) + self.assertEqual(rxn1.ts_species.ts_guesses[0].initial_xyz['symbols'], ('H', 'H', 'O')) + self.assert_h_abs_reactive_core(xyz=rxn1.ts_species.ts_guesses[0].initial_xyz, a=0, h=1, b=2, msg='H2 + O') expected_xyz = {'symbols': ('H', 'H', 'O'), 'isotopes': (1, 1, 16), - 'coords': ((0.0, 0.0, -1.8806939503689344), - (0.0, 0.0, -0.9839219710369642), (0.0, 0.0, 0.1804968455295033))} - self.assertTrue(almost_equal_coords(rxn1.ts_species.ts_guesses[0].initial_xyz, expected_xyz)) + 'coords': ((0.0, 0.0, -2.4256106790554215), + (0.0, 0.0, -1.1356106790554215), (0.0, 0.0, 0.22438932094457886))} + self.assertTrue(almost_equal_coords(rxn1.ts_species.ts_guesses[0].initial_xyz, expected_xyz, + rtol=1e-3, atol=1e-3)) # H + OH <=> H2 + O rxn2 = ARCReaction(r_species=[self.h, self.oh], p_species=[self.h2, self.o]) @@ -493,11 +601,7 @@ def test_heuristics_for_h_abstraction_1(self): heuristics_2.execute_incore() self.assertEqual(len(rxn2.ts_species.ts_guesses), 1) self.assertEqual(rxn2.ts_species.ts_guesses[0].initial_xyz['symbols'], ('H', 'O', 'H')) - expected_xyz = {'symbols': ('H', 'O', 'H'), 'isotopes': (1, 16, 1), - 'coords': ((0.0, 0.0, 1.8806939503689346), - (0.0, 0.0, -0.1804968455295031), - (0.0, 0.0, 0.9839219710369642))} - self.assertTrue(almost_equal_coords(rxn2.ts_species.ts_guesses[0].initial_xyz, expected_xyz)) + self.assert_h_abs_reactive_core(xyz=rxn2.ts_species.ts_guesses[0].initial_xyz, a=1, h=2, b=0, msg='H + OH') # OH + H <=> H2 + O rxn3 = ARCReaction(r_species=[self.oh, self.h], p_species=[self.h2, self.o]) @@ -511,11 +615,7 @@ def test_heuristics_for_h_abstraction_1(self): heuristics_3.execute_incore() self.assertEqual(len(rxn3.ts_species.ts_guesses), 1) self.assertEqual(rxn3.ts_species.ts_guesses[0].initial_xyz['symbols'], ('O', 'H', 'H')) - expected_xyz = {'symbols': ('O', 'H', 'H'), 'isotopes': (16, 1, 1), - 'coords': ((0.0, 0.0, -0.1804968455295031), - (0.0, 0.0, 0.9839219710369642), - (0.0, 0.0, 1.8806939503689346))} - self.assertTrue(almost_equal_coords(rxn3.ts_species.ts_guesses[0].initial_xyz, expected_xyz)) + self.assert_h_abs_reactive_core(xyz=rxn3.ts_species.ts_guesses[0].initial_xyz, a=0, h=1, b=2, msg='OH + H') # CH4 + H <=> CH3 + H2 ch4 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) @@ -540,14 +640,8 @@ def test_heuristics_for_h_abstraction_1(self): self.assertEqual(rxn4.ts_species.multiplicity, 2) self.assertEqual(len(rxn4.ts_species.ts_guesses), 4) # No dihedral scans for H attacking at 180 degrees. self.assertTrue(rxn4.ts_species.ts_guesses[0].success) - expected_xyz = {'symbols': ('C', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1, 1, 1), - 'coords': ((-0.14348351563387568, 3.0646463033967564e-08, 1.4446001062040636e-09), - (1.1671558180910826, -1.5372222794685086e-07, -1.0128514993379412e-07), - (-0.5075499920716767, -0.5148677050021889, -0.8917769786794892), - (-0.5075499920716767, -0.5148677050021889, 0.8917772176782378), - (-0.5075499920716767, 1.0297354786962867, 1.666119475718375e-08), - (2.063927797423041, -2.798718649055232e-07, -1.7157539600187732e-07))} - self.assertTrue(almost_equal_coords(rxn4.ts_species.ts_guesses[0].initial_xyz, expected_xyz)) + self.assertEqual(rxn4.ts_species.ts_guesses[0].initial_xyz['symbols'], ('C', 'H', 'H', 'H', 'H', 'H')) + self.assert_h_abs_reactive_core(xyz=rxn4.ts_species.ts_guesses[0].initial_xyz, a=0, h=1, b=5, msg='CH4 + H') def test_h_abstraction_with_empty_product_dicts(self): """ @@ -730,36 +824,26 @@ def test_heuristics_for_h_abstraction_5(self): self.assertEqual(rxn8.ts_species.charge, 0) self.assertEqual(rxn8.ts_species.multiplicity, 3) self.assertEqual(len(rxn8.ts_species.ts_guesses), 6) - self.assertTrue(almost_equal_coords(rxn8.ts_species.ts_guesses[0].initial_xyz, - {'symbols': ('N', 'C', 'O', 'N', 'H', 'H'), - 'isotopes': (14, 12, 16, 14, 1, 1), - 'coords': ( - (-3.657596721635545e-09, 0.08876698337705413, 0.9329034620293603), - (-3.657596721635545e-09, 0.8029731964901674, 0.005816759497990542), - (9.609228590831403e-09, 1.4947776654181317, -0.9420542025729876), - (-3.657596721635545e-09, -2.2452231277287726, 0.16101001965294726), - (-3.657596721635545e-09, -1.0762904363668742, 0.547597552628467), - (-3.657596721635545e-09, -2.2452231277287726, -0.8649899620469162))})) - self.assertTrue(almost_equal_coords(rxn8.ts_species.ts_guesses[1].initial_xyz, - {'symbols': ('N', 'C', 'O', 'N', 'H', 'H'), - 'isotopes': (14, 12, 16, 14, 1, 1), - 'coords': ((0.7304309896785263, 0.4813753237349452, -0.26044855417424406), - (-0.22605518455485318, 0.6753956743209286, 0.3853614572987496), - (-1.2013936879956266, 0.8535814918821567, 1.0130695121236126), - (0.7304309896785263, -1.8526147873708816, -1.032341996550657), - (0.7304309896785263, -0.6836820960089831, -0.6457544635751373), - (0.7304309896785263, -1.8526147873708816, - -2.0583419782505206))})) - self.assertTrue(almost_equal_coords(rxn8.ts_species.ts_guesses[2].initial_xyz, - {'symbols': ('N', 'C', 'O', 'N', 'H', 'H'), - 'isotopes': (14, 12, 16, 14, 1, 1), - 'coords': ( - (-0.7304309882994099, 0.48137532979234665, -0.26044855803662603), - (0.22605518593396912, 0.6753956803783296, 0.3853614534363681), - (1.201393684372416, 0.8535814759681686, 1.0130695222708512), - (-0.7304309882994099, -1.8526147813134801, -1.032342000413039), - (-0.7304309882994099, -0.6836820899515816, -0.6457544674375193), - (-0.7304309882994099, -1.8526147813134801, -2.0583419821129025))})) + seeds = h_abstraction(reaction=rxn8, dihedral_increment=120) + self.assertEqual(len(seeds), len(rxn8.ts_species.ts_guesses)) + for i, (seed, ts_guess) in enumerate(zip(seeds, rxn8.ts_species.ts_guesses)): + self.assertEqual(ts_guess.initial_xyz['symbols'], ('N', 'C', 'O', 'N', 'H', 'H')) + self.assertEqual(ts_guess.initial_xyz['isotopes'], (14, 12, 16, 14, 1, 1)) + self.assertTrue(almost_equal_coords(ts_guess.initial_xyz, seed['xyz'])) + reactive_atoms = seed['metadata']['reactive_atoms'] + self.assertEqual(reactive_atoms['A'], 3) + self.assertEqual(reactive_atoms['B'], 0) + self.assertIn(reactive_atoms['H'], [4, 5]) + self.assert_h_abs_reactive_core(xyz=ts_guess.initial_xyz, + a=reactive_atoms['A'], + h=reactive_atoms['H'], + b=reactive_atoms['B'], + msg=f'NCO + NH2 guess {i}', + ) + for i, j in itertools.combinations(range(len(rxn8.ts_species.ts_guesses)), 2): + self.assertFalse(almost_equal_coords(rxn8.ts_species.ts_guesses[i].initial_xyz, + rxn8.ts_species.ts_guesses[j].initial_xyz), + msg=f'NCO + NH2 guesses {i} and {j} are not unique') def test_heuristics_for_h_abstraction_6(self): # butenylnebzene + CCOO <=> butenylnebzene_rad + CCOOH @@ -949,6 +1033,19 @@ def test_heuristics_for_h_abstraction_8(self): h2o = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) rxn12 = ARCReaction(r_species=[nh3, oh], p_species=[nh2, h2o]) self.assertEqual(rxn12.family, 'H_Abstraction') + raw_seeds = h_abstraction(reaction=rxn12, dihedral_increment=60) + expected_reactive_atoms = [ + {'A': product_dict['r_label_map']['*1'], + 'H': product_dict['r_label_map']['*2'], + 'B': product_dict['r_label_map']['*3']} + for product_dict in rxn12.product_dicts + ] + for seed in raw_seeds: + reactive_atoms = seed['metadata']['reactive_atoms'] + self.assertIn(reactive_atoms, expected_reactive_atoms) + self.assertTrue(seed['xyz']['symbols'][reactive_atoms['H']].startswith('H')) + self.assertFalse(seed['xyz']['symbols'][reactive_atoms['A']].startswith('H')) + self.assertFalse(seed['xyz']['symbols'][reactive_atoms['B']].startswith('H')) heuristics_12 = HeuristicsAdapter(job_type='tsg', reactions=[rxn12], testing=True, @@ -1294,6 +1391,7 @@ def test_keeping_atom_order_in_ts(self): def test_combine_coordinates_with_redundant_atoms(self): """Test the combine_coordinates_with_redundant_atoms() function.""" + expected_symbols = ('C', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'C', 'H', 'H', 'H', 'H', 'H') ts_xyz = combine_coordinates_with_redundant_atoms( xyz_1=self.ccooh_xyz, xyz_2=self.c2h6_xyz, @@ -1309,6 +1407,27 @@ def test_combine_coordinates_with_redundant_atoms(self): d3=0, reactants_reversed=False, ) + self.assertEqual(ts_xyz['symbols'], expected_symbols) + self.assert_h_abs_reactive_core(xyz=ts_xyz, a=3, h=9, b=11, a2=180.0, msg='CCOOH + C2H6, a2 = 180') + + # An explicit multiplicative stretch reproduces the pre-absolute-distance geometry. + ts_xyz_stretched = combine_coordinates_with_redundant_atoms( + xyz_1=self.ccooh_xyz, + xyz_2=self.c2h6_xyz, + mol_1=ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz).mol, + mol_2=ARCSpecies(label='C2H6', smiles='CC', xyz=self.c2h6_xyz).mol, + reactant_2=ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz), + h1=9, + h2=5, + c=2, + d=0, + r1_stretch=1.2, + r2_stretch=1.2, + a2=180, + d2=None, + d3=0, + reactants_reversed=False, + ) expected_xyz = { 'symbols': ('C', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'C', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 16, 1, 1, 1, 1, 1, 1, 12, 12, 1, 1, 1, 1, 1), @@ -1329,7 +1448,7 @@ def test_combine_coordinates_with_redundant_atoms(self): (-2.5339216518555596, -0.010500763584261032, -0.3074446156924231), (-1.3956300152903744, -1.5298551170593468, 2.099824519093628), (-2.320130400140207, -0.017878639343882265, 2.183641197255159))} - self.assertTrue(almost_equal_coords(ts_xyz, expected_xyz)) + self.assertTrue(almost_equal_coords(ts_xyz_stretched, expected_xyz, rtol=1e-3, atol=1e-3)) ts_xyz = combine_coordinates_with_redundant_atoms(xyz_1=self.ccooh_xyz, xyz_2=self.c2h6_xyz, @@ -1348,6 +1467,29 @@ def test_combine_coordinates_with_redundant_atoms(self): d3=120, reactants_reversed=False, ) + self.assertEqual(ts_xyz['symbols'], expected_symbols) + self.assert_h_abs_reactive_core(xyz=ts_xyz, a=3, h=9, b=11, a2=150.0, msg='CCOOH + C2H6, a2 = 150') + + ts_xyz_stretched = combine_coordinates_with_redundant_atoms(xyz_1=self.ccooh_xyz, + xyz_2=self.c2h6_xyz, + mol_1=ARCSpecies(label='CCOOH', smiles='CCOO', + xyz=self.ccooh_xyz).mol, + mol_2=ARCSpecies(label='C2H6', smiles='CC', + xyz=self.c2h6_xyz).mol, + reactant_2=ARCSpecies(label='C2H5', + smiles='C[CH2]', + xyz=self.c2h5_xyz), + h1=9, + h2=5, + c=2, + d=0, + r1_stretch=1.2, + r2_stretch=1.2, + a2=150, + d2=30, + d3=120, + reactants_reversed=False, + ) expected_xyz = { 'symbols': ('C', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'C', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 16, 1, 1, 1, 1, 1, 1, 12, 12, 1, 1, 1, 1, 1), @@ -1368,7 +1510,7 @@ def test_combine_coordinates_with_redundant_atoms(self): (-0.6957522841048858, -2.2546439158856977, 1.860657306122405), (-2.5436577650309165, -0.1202574787360858, 0.6653623387328302), (-1.2481402862001663, -0.961958946530624, -0.20697484805157007))} - self.assertTrue(almost_equal_coords(ts_xyz, expected_xyz)) + self.assertTrue(almost_equal_coords(ts_xyz_stretched, expected_xyz, rtol=1e-3, atol=1e-3)) def test_get_new_zmat2_map(self): """Test the get_new_zmat_2_map() function.""" @@ -2255,6 +2397,69 @@ def test_check_ts_bonds(self): result = check_ts_bonds(initial_xyz, [7, 8, 9, 2, 4]) self.assertTrue(result) + def test_execute_incore_propagates_the_seed_method_and_family(self): + """A seed's ``method`` and ``family`` reach the TSGuess and the saved geometry filename.""" + project_directory = os.path.join(ARC_TESTING_PATH, 'heuristics_seed_plumbing') + self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True) + rxn = ARCReaction(r_species=[self.h2, self.o], p_species=[self.h, self.oh]) + seed_xyz = str_to_xyz("""H 0.0 0.0 -1.90 + H 0.0 0.0 -0.98 + O 0.0 0.0 0.18""") + seed = {'xyz': seed_xyz, + 'method': 'Heuristics-XY', + 'family': 'XY_Addition_MultipleBond', + 'source_adapter': 'heuristics', + 'metadata': {}} + heuristics = HeuristicsAdapter(job_type='tsg', + reactions=[rxn], + testing=True, + project='test', + project_directory=project_directory, + dihedral_increment=120, + ) + with patch('arc.job.adapters.ts.heuristics.get_ts_seeds', return_value=[seed]): + heuristics.execute_incore() + self.assertEqual(len(rxn.ts_species.ts_guesses), 1) + ts_guess = rxn.ts_species.ts_guesses[0] + self.assertEqual(ts_guess.method, 'heuristics-xy') + self.assertEqual(ts_guess.family, 'XY_Addition_MultipleBond') + self.assertEqual(ts_guess.method_index, 0) + self.assertTrue(ts_guess.success) + self.assertTrue(almost_equal_coords(ts_guess.initial_xyz, seed_xyz)) + self.assertTrue(os.path.isfile(os.path.join(heuristics.local_path, 'Heuristics-XY_0.xyz'))) + + def test_execute_incore_merges_a_duplicate_seed_into_an_existing_guess(self): + """A seed identical to an existing successful guess adds a source instead of a new guess.""" + project_directory = os.path.join(ARC_TESTING_PATH, 'heuristics_seed_merge') + self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True) + rxn = ARCReaction(r_species=[self.h2, self.o], p_species=[self.h, self.oh]) + rxn.ts_species = ARCSpecies(label='TS', is_ts=True, charge=rxn.charge, multiplicity=rxn.multiplicity) + seed_xyz = str_to_xyz("""H 0.0 0.0 -1.90 + H 0.0 0.0 -0.98 + O 0.0 0.0 0.18""") + existing = TSGuess(method='GCN', method_index=0, success=True, xyz=seed_xyz) + rxn.ts_species.append_ts_guess(existing) + seed = {'xyz': seed_xyz, + 'method': 'Heuristics', + 'family': rxn.family, + 'source_adapter': 'heuristics', + 'metadata': {}} + heuristics = HeuristicsAdapter(job_type='tsg', + reactions=[rxn], + testing=True, + project='test', + project_directory=project_directory, + dihedral_increment=120, + ) + with patch('arc.job.adapters.ts.heuristics.get_ts_seeds', return_value=[seed]): + heuristics.execute_incore() + self.assertEqual(len(rxn.ts_species.ts_guesses), 1) + merged = rxn.ts_species.ts_guesses[0] + self.assertEqual(merged.method, 'gcn') + self.assertIn('heuristics', merged.method_sources) + self.assertIn('gcn', merged.method_sources) + self.assertFalse(os.path.isfile(os.path.join(heuristics.local_path, 'Heuristics_0.xyz'))) + @classmethod def tearDownClass(cls): """ @@ -2265,5 +2470,332 @@ def tearDownClass(cls): shutil.rmtree(os.path.join(ARC_TESTING_PATH, sub), ignore_errors=True) +class TestHeuristicsHub(unittest.TestCase): + """Unit tests for shared heuristic seed and CREST-constraint helpers.""" + + def test_get_ts_seeds_h_abstraction(self): + rxn = SimpleNamespace(family='H_Abstraction') + with patch('arc.job.adapters.ts.heuristics.h_abstraction', + return_value=[{'xyz': {'symbols': ('H',), 'coords': ((0.0, 0.0, 0.0),), 'isotopes': (1,)}, + 'method': 'Heuristics', + 'metadata': {'source': 'reaction_mapping'}}]): + seeds = get_ts_seeds(reaction=rxn, base_adapter='heuristics', dihedral_increment=60) + self.assertEqual(len(seeds), 1) + self.assertEqual(seeds[0]['family'], 'H_Abstraction') + self.assertEqual(seeds[0]['method'], 'Heuristics') + self.assertEqual(seeds[0]['source_adapter'], 'heuristics') + self.assertEqual(seeds[0]['metadata'], {'source': 'reaction_mapping'}) + + def test_get_ts_seeds_hydrolysis(self): + rxn = SimpleNamespace(family='carbonyl_based_hydrolysis') + xyz = {'symbols': ('O',), 'coords': ((0.0, 0.0, 0.0),), 'isotopes': (16,)} + with patch('arc.job.adapters.ts.heuristics.hydrolysis', + return_value=([xyz], ['carbonyl_based_hydrolysis'], [[0, 1, 2]])): + seeds = get_ts_seeds(reaction=rxn, base_adapter='heuristics') + self.assertEqual(len(seeds), 1) + self.assertEqual(seeds[0]['family'], 'carbonyl_based_hydrolysis') + self.assertEqual(seeds[0]['xyz'], xyz) + self.assertEqual(seeds[0]['metadata'], {'indices': [0, 1, 2]}) + + def test_get_wrapper_constraints_crest(self): + rxn = SimpleNamespace(family='H_Abstraction') + xyz = str_to_xyz("""O 0.0000 0.0000 0.0000 + H 0.0000 0.0000 0.9600 + O 0.9000 0.0000 0.0000""") + seed = {'xyz': xyz, 'family': rxn.family} + constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) + self.assertIsInstance(constraints, dict) + self.assertLessEqual({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'}, set(constraints)) + self.assertEqual( + (constraints['A'], constraints['H'], constraints['B']), + constraints['angle_atoms'], + ) + self.assertEqual(len(constraints['atoms']), 3) + self.assertEqual(len(constraints['distance_pairs']), 2) + + def test_get_wrapper_constraints_crest_symmetric_oh_oh(self): + """The transferring H must be bracketed by the two O atoms in either atom ordering.""" + rxn = SimpleNamespace(family='H_Abstraction') + seeds_and_expected_atoms = [ + (str_to_xyz("""O 0.00000000 -0.02752832 -1.20590500 + H 0.00000000 -0.02752832 -0.03383145 + O 0.00000000 -0.02752832 1.12142787 + H 0.00000000 0.90131726 1.37454478"""), 1), + (str_to_xyz("""O 0.00000000 -0.02752832 1.12142787 + H 0.00000000 0.90131726 1.37454478 + O 0.00000000 -0.02752832 -1.20590500 + H 0.00000000 -0.02752832 -0.03383145"""), 3), + ] + for xyz, expected_h_atom in seeds_and_expected_atoms: + with self.subTest(symbols=xyz['symbols']): + seed = { + 'xyz': xyz, + 'family': rxn.family, + 'metadata': {'reactive_atoms': {'A': 0, 'H': expected_h_atom, 'B': 2}}, + } + constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) + self.assertEqual(constraints['angle_atoms'][1], expected_h_atom) + self.assertSetEqual({constraints['angle_atoms'][0], constraints['angle_atoms'][2]}, {0, 2}) + self.assertFalse(xyz['symbols'][constraints['angle_atoms'][0]].startswith('H')) + self.assertFalse(xyz['symbols'][constraints['angle_atoms'][2]].startswith('H')) + + def test_get_wrapper_constraints_crest_supports_abstraction_by_an_h_atom(self): + """R-H + H <=> R + H2: the acceptor may itself be a hydrogen. + + Only the transferred atom has to be a hydrogen. CREST constraints are index-based, so + requiring a heavy acceptor silently excluded every abstraction by an H atom. + """ + rxn = SimpleNamespace(family='H_Abstraction') + xyz = str_to_xyz("""C 0.0000 0.0000 0.0000 + H 1.0300 0.0000 -0.3600 + H -0.5100 -0.8900 -0.3600 + H -0.5100 0.8900 -0.3600 + H 0.0000 0.0000 1.3000 + H 0.0000 0.0000 2.2000""") + seed = {'xyz': xyz, 'family': rxn.family, + 'metadata': {'reactive_atoms': {'A': 0, 'H': 4, 'B': 5}}} + constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) + self.assertIsNotNone(constraints) + self.assertEqual(constraints['atoms'], (0, 4, 5)) + self.assertEqual(constraints['distance_pairs'], ((0, 4), (4, 5))) + self.assertEqual(constraints['angle_atoms'], (0, 4, 5)) + + def test_get_wrapper_constraints_crest_rejects_hydrogen_as_heavy_atom(self): + rxn = SimpleNamespace(family='H_Abstraction') + xyz = str_to_xyz("""O -1.0000 0.0000 0.0000 + H 0.0000 0.0000 0.0000 + O 1.0000 0.0000 0.0000 + H 2.0000 0.0000 0.0000""") + seed = { + 'xyz': xyz, + 'family': rxn.family, + 'metadata': {'reactive_atoms': {'A': 0, 'H': 1, 'B': 3}}, + } + self.assertIsNone(get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)) + + def test_get_ts_seeds_preserves_invalid_explicit_reactive_atoms(self): + """Do not hide an invalid generator mapping with the geometric compatibility fallback.""" + rxn = SimpleNamespace(family='H_Abstraction') + xyz = str_to_xyz("""O -1.0000 0.0000 0.0000 + H 0.0000 0.0000 0.0000 + O 1.0000 0.0000 0.0000 + H 2.0000 0.0000 0.0000""") + invalid_atoms = {'A': 0, 'H': 1, 'B': 3} + with patch('arc.job.adapters.ts.heuristics.h_abstraction', return_value=[{ + 'xyz': xyz, + 'method': 'Heuristics', + 'metadata': {'reactive_atoms': invalid_atoms}, + }]): + seed = get_ts_seeds(reaction=rxn, base_adapter='heuristics')[0] + self.assertEqual(seed['metadata']['reactive_atoms'], invalid_atoms) + self.assertIsNone(get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)) + + def test_get_wrapper_constraints_crest_prefers_explicit_reactive_atoms(self): + """Explicit generator metadata wins when distances would select a spectator hydrogen.""" + rxn = SimpleNamespace(family='H_Abstraction') + xyz = str_to_xyz("""O -1.0000 0.0000 0.0000 + H 0.0000 0.0000 0.0000 + O 1.0000 0.0000 0.0000 + H 0.0000 1.5000 0.0000""") + seed = { + 'xyz': xyz, + 'family': rxn.family, + 'metadata': {'reactive_atoms': {'A': 0, 'H': 3, 'B': 2}}, + } + self.assertEqual( + get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed), + { + 'A': 0, + 'H': 3, + 'B': 2, + 'atoms': (0, 3, 2), + 'distance_pairs': ((0, 3), (3, 2)), + 'angle_atoms': (0, 3, 2), + }, + ) + + def test_get_wrapper_constraints_crest_xy_addition(self): + """XY constraints follow the exact family-label mapping for both seed orderings.""" + rxn = SimpleNamespace(family='XY_Addition_MultipleBond') + xyz = str_to_xyz("""C 0.0 0.0 0.667 + C 0.0 0.0 -0.667 + H 0.0 1.6 0.667 + Cl 0.0 2.1 -0.667""") + for reactive_atoms in ( + {'*1': 1, '*2': 0, '*3': 2, '*4': 3}, + {'*1': 0, '*2': 1, '*3': 2, '*4': 3}, + ): + with self.subTest(reactive_atoms=reactive_atoms): + seed = { + 'xyz': xyz, + 'family': rxn.family, + 'metadata': {'reactive_atoms': reactive_atoms}, + } + self.assertEqual( + get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed), + { + 'atoms': tuple(reactive_atoms[label] for label in ('*1', '*2', '*3', '*4')), + 'distance_pairs': ( + (reactive_atoms['*1'], reactive_atoms['*3']), + (reactive_atoms['*2'], reactive_atoms['*4']), + (reactive_atoms['*3'], reactive_atoms['*4']), + ), + }, + ) + + def test_get_wrapper_constraints_crest_rejects_invalid_explicit_xy_atoms(self): + """Do not infer an XY mapping when explicit generator metadata is invalid.""" + rxn = SimpleNamespace(family='XY_Addition_MultipleBond') + xyz = str_to_xyz("""C 0.0 0.0 0.667 + C 0.0 0.0 -0.667 + H 0.0 1.6 0.667 + Cl 0.0 2.1 -0.667""") + invalid_atoms = {'*1': 0, '*2': 1, '*3': 2, '*4': 2} + with patch('arc.job.adapters.ts.xy_addition.xy_addition', return_value=[{ + 'xyz': xyz, + 'method': 'Heuristics-XY', + 'metadata': {'reactive_atoms': invalid_atoms}, + }]): + seed = get_ts_seeds(reaction=rxn)[0] + self.assertEqual(seed['metadata']['reactive_atoms'], invalid_atoms) + self.assertIsNone(get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)) + + def test_get_wrapper_constraints_crest_unsupported_family(self): + rxn = SimpleNamespace(family='carbonyl_based_hydrolysis') + xyz = str_to_xyz("""O 0.0000 0.0000 0.0000 + H 0.0000 0.0000 0.9600 + H 0.9000 0.0000 0.0000""") + seed = {'xyz': xyz, 'family': rxn.family} + atoms = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) + self.assertIsNone(atoms) + + def test_h_atom_bond_cutoff_brackets_a_bound_hydrogen(self): + """An H is a valid free H just past a carbon's bonding cutoff and bound just below it. + + Carbon's cutoff is ``H_ATOM_BOND_CUTOFF``, the floor of every element-specific cutoff, and + the two cases sit 0.05 A on either side of it, so widening or narrowing it changes one of + them. + """ + self.assertAlmostEqual(H_ATOM_BOND_CUTOFF, 1.5, delta=1e-9) + self.assertAlmostEqual(_h_atom_bond_cutoff('C'), H_ATOM_BOND_CUTOFF, delta=1e-9) + for neighbour_distance, expected in ((H_ATOM_BOND_CUTOFF - 0.05, False), + (H_ATOM_BOND_CUTOFF + 0.05, True)): + with self.subTest(neighbour_distance=neighbour_distance): + xyz = str_to_xyz(f"""H 0.0 0.0 0.0 + H 0.0 0.0 1.6 + O 0.0 0.0 2.6 + C {neighbour_distance} 0.0 0.0""") + self.assertEqual(_is_valid_h_abs_atom_assignment(xyz=xyz, atoms={'A': 0, 'H': 1, 'B': 2}), + expected) + + def test_h_atom_bond_cutoff_is_element_aware(self): + """The bonding cutoff widens with the estimated X--H bond length, never below the floor. + + The rows cover each source the cutoff can draw on: the floor, ARC's tabulated single bond + length, the covalent radii sum, and the largest radius of an element that is tabulated only + per spin state. + """ + for symbol, expected in (('H', 1.5), + ('C', 1.5), + ('Cl', 1.5), + ('Si', 1.702), + ('I', 1.8745), + ('Ge', 1.7365), + ('Sn', 1.955), + ('Fe', 2.1045), + ): + with self.subTest(symbol=symbol): + self.assertAlmostEqual(_h_atom_bond_cutoff(symbol), expected, places=4) + for symbol in ('Xx', 26): + with self.subTest(symbol=symbol): + self.assertAlmostEqual(_h_atom_bond_cutoff(symbol), H_ATOM_BOND_CUTOFF, delta=1e-9) + + def test_a_long_heavy_element_hydride_bond_is_not_a_free_h(self): + """A hydrogen at a Ge--H, Sn--H, I--H or Fe--H bond length is bound, not a free H radical. + + All four bonds are longer than ``H_ATOM_BOND_CUTOFF``, so a single flat cutoff reads them + as free hydrogens. The distances are experimental bond lengths, and the free counter-case + sits at a non-bonded separation from the same element. + """ + for symbol, bond_length, free_distance in (('Ge', 1.53, 2.50), + ('Sn', 1.70, 2.70), + ('I', 1.61, 2.60), + ('Fe', 1.55, 2.60), + ): + for neighbour_distance, expected in ((bond_length, False), (free_distance, True)): + with self.subTest(symbol=symbol, neighbour_distance=neighbour_distance): + xyz = str_to_xyz(f"""H 0.0 0.0 0.0 + H 0.0 0.0 1.6 + O 0.0 0.0 2.6 + {symbol} {neighbour_distance} 0.0 0.0""") + self.assertEqual(_is_valid_h_abs_atom_assignment(xyz=xyz, atoms={'A': 0, 'H': 1, 'B': 2}), + expected) + + def test_h_abs_atom_selection_tie_break_keeps_the_first_hydrogen(self): + """When two hydrogens tie on the sum of the distances to their partners, the first is kept. + + Each hydrogen is a free H atom, so the candidate partners of either are the other hydrogen + and the two oxygens. Both hydrogens reach the same partner-distance sum of 4.236 A, and the + lower index wins. + """ + xyz = str_to_xyz("""O -2.0 0.0 0.0 + O 2.0 0.0 0.0 + H 0.0 1.0 0.0 + H 0.0 -1.0 0.0""") + self.assertEqual(_get_h_abs_atoms_from_xyz(xyz), {'H': 2, 'A': 3, 'B': 0}) + + def test_h_abs_atoms_resolve_a_free_h_reactant(self): + """A free H atom is a valid H-abstraction reactant, as in CH4 + H <=> CH3 + H2. + + The only heavy atom of this TS is the methyl carbon, so searching for the two partners of + the transferred hydrogen among the heavy atoms alone leaves fewer than two candidates and + resolves no triad at all. The three methyl hydrogens are bound to the carbon and stay + excluded. + """ + xyz = str_to_xyz("""C 0.000 0.000 0.000 + H 0.000 0.000 1.290 + H 1.028 0.000 -0.363 + H -0.514 0.890 -0.363 + H -0.514 -0.890 -0.363 + H 0.000 0.000 2.650""") + self.assertEqual(_get_h_abs_atoms_from_xyz(xyz), {'H': 1, 'A': 0, 'B': 5}) + self.assertTrue(_is_valid_h_abs_atom_assignment(xyz=xyz, atoms={'A': 0, 'H': 1, 'B': 5})) + rxn = SimpleNamespace(family='H_Abstraction') + constraints = get_wrapper_constraints(wrapper='crest', + reaction=rxn, + seed={'xyz': xyz, 'family': 'H_Abstraction'}, + ) + self.assertEqual(constraints['atoms'], (0, 1, 5)) + self.assertEqual(constraints['distance_pairs'], ((0, 1), (1, 5))) + self.assertEqual(constraints['angle_atoms'], (0, 1, 5)) + + def test_get_h_abs_atoms_from_xyz_without_two_partner_atoms(self): + """A geometry with fewer than two candidate partners has no A--H--B triad. + + Water has three atoms and still yields no triad: each of its hydrogens is bound to the + oxygen, so it is a spectator rather than a free H atom that the other could take as a + partner. The two hydrogens of this equilibrium geometry are 1.514 A apart, past the cutoff + of a hydrogen neighbour, so it is the oxygen at 0.958 A that binds them. + """ + xyz = str_to_xyz("""O 0.0 0.0 0.0 + H 0.0 0.0 0.97""") + self.assertIsNone(_get_h_abs_atoms_from_xyz(xyz)) + water = str_to_xyz("""O 0.0000 0.0000 0.1173 + H 0.0000 0.7572 -0.4692 + H 0.0000 -0.7572 -0.4692""") + self.assertGreater(calculate_distance(coords=water['coords'], atoms=[1, 2]), H_ATOM_BOND_CUTOFF) + self.assertIsNone(_get_h_abs_atoms_from_xyz(water)) + + def test_get_ts_seeds_unsupported_adapter(self): + rxn = SimpleNamespace(family='H_Abstraction') + with self.assertRaises(ValueError): + get_ts_seeds(reaction=rxn, base_adapter='gcn') + + def test_get_wrapper_constraints_unsupported_wrapper(self): + rxn = SimpleNamespace(family='H_Abstraction') + with self.assertRaises(ValueError): + get_wrapper_constraints(wrapper='foo_wrapper', reaction=rxn, seed={}) + + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/arc/job/adapters/ts/seed_hub.py b/arc/job/adapters/ts/seed_hub.py new file mode 100644 index 0000000000..4191702fd8 --- /dev/null +++ b/arc/job/adapters/ts/seed_hub.py @@ -0,0 +1,505 @@ +""" +Shared TS-seed and wrapper-constraint hub. + +This module centralizes: +1. How TS seeds are requested from a base TS-search adapter. +2. How wrapper adapters (e.g., CREST) request family-specific constraints for a seed. +""" + +from functools import lru_cache +from typing import Dict, List, Optional, Sequence + +from arc.common import (COVALENT_RADII, + SINGLE_BOND_LENGTH, + almost_equal_coords, + get_atom_radius, + get_logger, + get_single_bond_length, + ) +from arc.species.converter import xyz_to_dmat + +logger = get_logger() + +# Shortest distance at which a hydrogen is taken to be free of a neighbour, used to tell a free H +# radical (a valid reactant) from a bound spectator hydrogen. Serves as the floor of every +# element-specific cutoff, and as the cutoff of an element whose bond to hydrogen cannot be +# estimated at all. The floor covers the H, B, C, N, O and F cutoffs, whose bonds to hydrogen are +# short enough that the slack alone would reject a TS-stretched one. +H_ATOM_BOND_CUTOFF = 1.5 + +# Factor applied to the estimated single bond length of an element with hydrogen to obtain that +# element's hydrogen bonding cutoff, where that exceeds the floor. It admits the stretching of a +# bond in a TS geometry. +H_ATOM_BOND_CUTOFF_SLACK = 1.15 + +# The reaction families handled by the heuristics hydrolysis seed builder, mirroring +# ``heuristics.FAMILY_SETS``. Held here so that resolving them does not import heuristics. +HYDROLYSIS_FAMILIES = ('carbonyl_based_hydrolysis', + 'ether_hydrolysis', + 'nitrile_hydrolysis', + ) + +# The hydrolysis mechanism roles whose mutual geometry CREST has to hold fixed: ``a`` is the +# electrophilic centre, ``b`` the leaving group, ``o`` the water oxygen and ``h1`` the water +# hydrogen that transfers from ``o`` to ``b``. +HYDROLYSIS_REACTIVE_ROLES = ('a', 'b', 'o', 'h1') + +# The order in which ``heuristics.hydrolysis()`` emits its per-guess index sequence. +HYDROLYSIS_INDEX_ORDER = ('a', 'b', 'e', 'o', 'h1', 'd') + +# All six pairwise separations among the four reactive atoms. Four atoms have six internal +# degrees of freedom, so the six pairwise distances determine the reactive core rigidly, up to +# a reflection that the reference geometry resolves. +HYDROLYSIS_DISTANCE_ROLE_PAIRS = (('a', 'b'), + ('a', 'o'), + ('a', 'h1'), + ('b', 'o'), + ('b', 'h1'), + ('o', 'h1'), + ) + + +def get_ts_seeds(reaction: 'ARCReaction', + base_adapter: str = 'heuristics', + dihedral_increment: Optional[int] = None, + ) -> List[dict]: + """ + Return TS seed entries from a base TS-search adapter. + + Seed schema: + - ``xyz`` (dict): Cartesian coordinates. + - ``family`` (str): The family associated with this seed. + - ``method`` (str): Human-readable generator label. + - ``source_adapter`` (str): Adapter id that generated the seed. + - ``metadata`` (dict, optional): Adapter-specific auxiliary fields. + + Args: + reaction: The ARC reaction object. + base_adapter: The underlying TS-search adapter providing seeds. + dihedral_increment: Optional scan increment used by adapters that support it. + """ + adapter = (base_adapter or '').lower() + if adapter != 'heuristics': + raise ValueError(f'Unsupported TS seed base adapter: {base_adapter}') + + # Lazily import to avoid circular imports with heuristics.py. + from arc.job.adapters.ts.heuristics import h_abstraction, hydrolysis + + xyz_entries = list() + if reaction.family == 'H_Abstraction': + xyzs = h_abstraction(reaction=reaction, dihedral_increment=dihedral_increment) + for entry in xyzs: + xyz = entry.get('xyz') if isinstance(entry, dict) else entry + method = entry.get('method', 'Heuristics') if isinstance(entry, dict) else 'Heuristics' + if xyz is not None: + entry_metadata = entry.get('metadata') if isinstance(entry, dict) else None + metadata = entry_metadata.copy() if isinstance(entry_metadata, dict) else {} + if 'reactive_atoms' not in metadata: + reactive_atoms = _get_h_abs_atoms_from_xyz(xyz) + if reactive_atoms is not None: + metadata['reactive_atoms'] = reactive_atoms + xyz_entries.append({ + 'xyz': xyz, + 'method': method, + 'family': reaction.family, + 'source_adapter': 'heuristics', + 'metadata': metadata, + }) + elif reaction.family in get_hydrolysis_families(): + try: + xyzs_raw, families, indices = hydrolysis(reaction=reaction) + except ValueError: + xyz_entries = list() + else: + for xyz, family, idx in zip(xyzs_raw, families, indices): + metadata = {'indices': idx} + reactive_atoms = get_hydrolysis_reactive_atoms(idx) + if reactive_atoms is not None: + metadata['reactive_atoms'] = reactive_atoms + xyz_entries.append({ + 'xyz': xyz, + 'method': 'Heuristics', + 'family': family, + 'source_adapter': 'heuristics', + 'metadata': metadata, + }) + elif reaction.family == 'XY_Addition_MultipleBond': + # Lazily import to keep the family-specific builder decoupled from this hub. + from arc.job.adapters.ts.xy_addition import xy_addition + for entry in xy_addition(reaction=reaction): + xyz_entries.append({ + 'xyz': entry['xyz'], + 'method': entry.get('method', 'Heuristics-XY'), + 'family': reaction.family, + 'source_adapter': 'heuristics', + 'metadata': entry.get('metadata', {}).copy(), + }) + return xyz_entries + + +def get_hydrolysis_families() -> List[str]: + """ + Return the reaction family names handled by the heuristics hydrolysis seed builder. + + These mirror ``heuristics.FAMILY_SETS`` and are held here so that resolving them does not + import ``arc.job.adapters.ts.heuristics``. ``test_hydrolysis_families_match_family_sets`` + fails if the two ever disagree. + + Returns: + List[str]: The families of both hydrolysis parameter sets. + """ + return list(HYDROLYSIS_FAMILIES) + + +def get_hydrolysis_reactive_atoms(indices) -> Optional[Dict[str, int]]: + """ + Return the named hydrolysis reactive atoms of a seed's raw index metadata. + + Two shapes are accepted: a mapping already keyed by the mechanism roles, of which only + ``a``, ``b``, ``o`` and ``h1`` are kept, and the positional sequence + ``(a, b, e, o, h1, d)`` that :func:`arc.job.adapters.ts.heuristics.hydrolysis` emits. + + Args: + indices: The ``indices`` entry of a hydrolysis seed's metadata. + + Returns: + Optional[Dict[str, int]]: The ``a``, ``b``, ``o`` and ``h1`` atom indices, or ``None`` + if the roles cannot be resolved from ``indices``. + """ + positional = None + if isinstance(indices, dict): + positional = indices + elif isinstance(indices, (list, tuple)) and len(indices) == len(HYDROLYSIS_INDEX_ORDER): + positional = dict(zip(HYDROLYSIS_INDEX_ORDER, indices)) + if positional is None or any(role not in positional for role in HYDROLYSIS_REACTIVE_ROLES): + return None + return {role: positional[role] for role in HYDROLYSIS_REACTIVE_ROLES} + + +def get_backup_ts_seeds(reaction: 'ARCReaction', + exclude_method: str = 'crest', + ) -> List[dict]: + """ + Build CREST seed entries from TS guesses that OTHER adapters already produced. + + This is a fallback for when CREST's own heuristic seed construction + (:func:`get_ts_seeds`) yields nothing -- e.g. a linear/cumulene reactive center + (such as HCCO in H_Abstraction) that the heuristic Z-matrix builder cannot + assemble. Any successful non-CREST TS guess already present on + ``reaction.ts_species.ts_guesses`` is a valid CREST seed: CREST only needs a seed + geometry plus the family reactive-atom constraints, and the constraints are + re-derived from the seed geometry by :func:`get_wrapper_constraints` -- they do not + depend on how the seed geometry was originally built. + + Seeds are returned with empty ``metadata`` so the wrapper-constraint derivation + re-infers the reactive atoms from the geometry itself (robust to the source + adapter's atom ordering). Guesses whose method contains ``exclude_method`` are + skipped so CREST is never seeded from a prior CREST result (feedback-loop guard). + + Args: + reaction: The ARC reaction object. + exclude_method: A method substring to exclude (default ``'crest'``). + + Returns: + List[dict]: Seed entries in the same schema as :func:`get_ts_seeds`. + """ + ts_species = getattr(reaction, 'ts_species', None) + ts_guesses = getattr(ts_species, 'ts_guesses', None) or list() + exclude = (exclude_method or '').lower() + seeds = list() + seen_xyzs = list() + for tsg in ts_guesses: + method = (getattr(tsg, 'method', '') or '').lower() + if not getattr(tsg, 'success', False): + continue + if exclude and exclude in method: + continue + xyz = getattr(tsg, 'opt_xyz', None) or getattr(tsg, 'initial_xyz', None) + if not isinstance(xyz, dict) or not xyz.get('symbols'): + continue + if any(almost_equal_coords(xyz, seen) for seen in seen_xyzs): + continue + seen_xyzs.append(xyz) + seeds.append({ + 'xyz': xyz, + 'method': getattr(tsg, 'method', None) or 'external', + 'family': reaction.family, + 'source_adapter': method or 'external', + 'metadata': {}, + }) + return seeds + + +def get_wrapper_constraints(wrapper: str, + reaction: 'ARCReaction', + seed: dict, + ) -> Optional[dict]: + """ + Return wrapper-specific constraints for a TS seed. + + Args: + wrapper: Wrapper adapter id (e.g., ``crest``). + reaction: The ARC reaction object. + seed: A seed entry returned by :func:`get_ts_seeds`. + """ + wrapper_name = (wrapper or '').lower() + if wrapper_name != 'crest': + raise ValueError(f'Unsupported wrapper adapter: {wrapper}') + return _get_crest_constraints(reaction=reaction, seed=seed) + + +def _get_crest_constraints(reaction: 'ARCReaction', seed: dict) -> Optional[dict]: + """ + Return a generic CREST constraint specification for a seed. + + The specification contains zero-based participating ``atoms`` and ``distance_pairs``. + H-abstraction additionally supplies ``angle_atoms`` so completed geometries retain + the seed's heavy-atom--H--heavy-atom orientation. + + The hydrolysis families pin all six pairwise distances among the four reactive atoms + ``a``, ``b``, ``o`` and ``h1``, which are read from the seed metadata: either from an + explicit ``reactive_atoms`` mapping or, failing that, from the raw ``indices`` entry. + ``None`` is returned when neither resolves to four distinct, in-range atoms whose ``h1`` + is a hydrogen. + """ + family = seed.get('family') or reaction.family + xyz = seed.get('xyz') + if xyz is None: + return None + metadata = seed.get('metadata') + explicit_atoms = metadata.get('reactive_atoms') if isinstance(metadata, dict) else None + if family == 'H_Abstraction': + reactive_atoms = explicit_atoms if explicit_atoms is not None else _get_h_abs_atoms_from_xyz(xyz) + if _is_valid_h_abs_atom_assignment(xyz=xyz, atoms=reactive_atoms): + return { + 'A': reactive_atoms['A'], + 'H': reactive_atoms['H'], + 'B': reactive_atoms['B'], + 'atoms': tuple(reactive_atoms[key] for key in ('A', 'H', 'B')), + 'distance_pairs': ( + (reactive_atoms['A'], reactive_atoms['H']), + (reactive_atoms['H'], reactive_atoms['B']), + ), + 'angle_atoms': tuple(reactive_atoms[key] for key in ('A', 'H', 'B')), + } + if explicit_atoms is not None: + logger.warning(f'Invalid explicit CREST H-abstraction atom assignment: {explicit_atoms}') + return None + if family == 'XY_Addition_MultipleBond': + if _is_valid_xy_atom_assignment(xyz=xyz, atoms=explicit_atoms): + return { + 'atoms': tuple(explicit_atoms[label] for label in ('*1', '*2', '*3', '*4')), + 'distance_pairs': ( + (explicit_atoms['*1'], explicit_atoms['*3']), + (explicit_atoms['*2'], explicit_atoms['*4']), + (explicit_atoms['*3'], explicit_atoms['*4']), + ), + } + logger.warning(f'Invalid explicit CREST XY-addition atom assignment: {explicit_atoms}') + return None + if family in get_hydrolysis_families(): + raw_indices = metadata.get('indices') if isinstance(metadata, dict) else None + reactive_atoms = explicit_atoms if explicit_atoms is not None \ + else get_hydrolysis_reactive_atoms(raw_indices) + if _is_valid_hydrolysis_atom_assignment(xyz=xyz, atoms=reactive_atoms): + return { + 'atoms': tuple(reactive_atoms[role] for role in HYDROLYSIS_REACTIVE_ROLES), + 'distance_pairs': tuple((reactive_atoms[role_1], reactive_atoms[role_2]) + for role_1, role_2 in HYDROLYSIS_DISTANCE_ROLE_PAIRS), + } + logger.warning(f'Could not determine the CREST hydrolysis reactive atoms of a {family} seed ' + f'from the metadata {metadata!r}, skipping this seed.') + return None + + +def _is_valid_hydrolysis_atom_assignment(xyz: dict, atoms: Optional[Dict[str, int]]) -> bool: + """ + Return whether ``atoms`` identifies four distinct, in-range hydrolysis reactive atoms. + + The transferred atom ``h1`` has to be a hydrogen, which is what distinguishes a correctly + ordered role assignment from one whose roles were permuted. + + Args: + xyz (dict): The seed geometry. + atoms (Optional[Dict[str, int]]): The candidate role-to-index mapping. + + Returns: + bool: Whether ``atoms`` is a usable hydrolysis reactive-atom assignment. + """ + symbols = xyz.get('symbols') if isinstance(xyz, dict) else None + if not symbols or not isinstance(atoms, dict) or set(atoms) != set(HYDROLYSIS_REACTIVE_ROLES): + return False + indices = tuple(atoms[role] for role in HYDROLYSIS_REACTIVE_ROLES) + if not all(isinstance(index, int) and not isinstance(index, bool) and 0 <= index < len(symbols) + for index in indices): + return False + if len(set(indices)) != len(HYDROLYSIS_REACTIVE_ROLES): + return False + return symbols[atoms['h1']] == 'H' + + +def _is_valid_xy_atom_assignment(xyz: dict, atoms: Optional[Dict[str, int]]) -> bool: + """Return whether ``atoms`` identifies four distinct, in-range XY recipe atoms.""" + symbols = xyz.get('symbols') if isinstance(xyz, dict) else None + if not symbols or not isinstance(atoms, dict) or set(atoms) != {'*1', '*2', '*3', '*4'}: + return False + indices = tuple(atoms[label] for label in ('*1', '*2', '*3', '*4')) + return (all(isinstance(index, int) and 0 <= index < len(symbols) for index in indices) + and len(set(indices)) == 4) + + +@lru_cache(maxsize=None) +def _h_atom_bond_cutoff(symbol: str) -> float: + """ + Return the distance below which ``symbol`` is taken to be covalently bound to a hydrogen. + + The cutoff is ``H_ATOM_BOND_CUTOFF_SLACK`` times an estimate of the ``symbol``--hydrogen single + bond length, floored at ``H_ATOM_BOND_CUTOFF``. The estimate is ARC's tabulated single bond + length where the pair has one, and the covalent radii sum otherwise, falling back to the largest + tabulated radius of an element that ``COVALENT_RADII`` keys only by hybridisation or spin state. + ``H_ATOM_BOND_CUTOFF`` is returned for a ``symbol`` with neither. + + Args: + symbol (str): The element symbol of the neighbour. + + Returns: + float: The bonding cutoff in Angstrom. + """ + if not isinstance(symbol, str): + return H_ATOM_BOND_CUTOFF + if f'{symbol}_H' in SINGLE_BOND_LENGTH or f'H_{symbol}' in SINGLE_BOND_LENGTH: + return max(H_ATOM_BOND_CUTOFF, H_ATOM_BOND_CUTOFF_SLACK * get_single_bond_length(symbol, 'H')) + radius = get_atom_radius(symbol) + if radius is None: + radius = max((value for key, value in COVALENT_RADII.items() if key.split('_')[0] == symbol), + default=None) + if radius is None: + return H_ATOM_BOND_CUTOFF + return max(H_ATOM_BOND_CUTOFF, H_ATOM_BOND_CUTOFF_SLACK * (radius + get_atom_radius('H'))) + + +def _is_free_h_atom(symbols: Sequence[str], + dmat, + index: int, + transferred_h_index: int, + ) -> bool: + """ + Return whether atom ``index`` is a hydrogen that is not covalently bound to any atom other than + the transferred hydrogen ``transferred_h_index``. + + A free H atom is a valid H-abstraction reactant, as in ``R-H + H <=> R + H2``, while a hydrogen + that is bound to something else is a spectator. Each neighbour is compared against the cutoff of + its own element, :func:`_h_atom_bond_cutoff`, so that a long bond such as Sn--H or I--H is + recognised as one. + + The neighbour scan skips both ``index`` itself and ``transferred_h_index``, and nothing else: an + ``index`` whose only close contact is the hypothesised transferred hydrogen is free, while every + other contact, including a TS-stretched one, keeps it bound. + + Args: + symbols (Sequence[str]): The element symbols of the geometry. + dmat: The interatomic distance matrix of the geometry. + index (int): The index of the atom to check. + transferred_h_index (int): The index of the transferred hydrogen. + + Returns: + bool: Whether ``index`` is a free hydrogen atom. + """ + if symbols[index] != 'H': + return False + return all(dmat[index][i] >= _h_atom_bond_cutoff(symbols[i]) + for i in range(len(symbols)) if i not in (index, transferred_h_index)) + + +def _is_valid_h_abs_atom_assignment(xyz: dict, atoms: Optional[Dict[str, int]]) -> bool: + """ + Return whether ``atoms`` identifies an ``A``--``H``--``B`` triad in ``xyz``, where ``H`` is the + transferred hydrogen and ``A`` and ``B`` are its two partner atoms. + + Only ``H`` has to be a hydrogen. ``A`` and ``B`` may be hydrogens themselves, as in + ``R-H + H <=> R + H2``: the CREST constraints are index-based, so the element of ``A`` and ``B`` + is irrelevant to them. Requiring both to be heavy atoms silently excluded every abstraction by + an H atom. + """ + symbols = xyz.get('symbols') if isinstance(xyz, dict) else None + if not symbols or not isinstance(atoms, dict) or set(atoms) != {'A', 'H', 'B'}: + return False + if any(not isinstance(atoms[key], int) or not 0 <= atoms[key] < len(symbols) for key in atoms): + return False + if len({atoms['A'], atoms['H'], atoms['B']}) != 3: + return False + if symbols[atoms['H']] != 'H': + return False + if all(symbols[atoms[key]] != 'H' for key in ('A', 'B')): + return True + # A hydrogen ``A`` or ``B`` is only meaningful when it is a free H atom, as in + # ``R-H + H <=> R + H2``. A hydrogen that is covalently bound to something else is a spectator + # that a faulty generator mapping named by mistake, which is what the heavy-atom-only rule was + # really guarding against. + dmat = xyz_to_dmat(xyz) + if dmat is None: + return False + for key in ('A', 'B'): + index = atoms[key] + if symbols[index] == 'H' and not _is_free_h_atom(symbols=symbols, + dmat=dmat, + index=index, + transferred_h_index=atoms['H'], + ): + return False + return True + + +def _get_h_abs_atoms_from_xyz(xyz: dict) -> Optional[Dict[str, int]]: + """ + Determine H-abstraction atoms from a TS guess. + + ``H`` is the transferred hydrogen, ``A`` the atom nearest to it and ``B`` the second nearest, + both chosen among the heavy atoms and the free hydrogens of the geometry. Free hydrogens are + identified with :func:`_is_free_h_atom`, the same check :func:`_is_valid_h_abs_atom_assignment` + applies, so that an abstraction by or from a free H atom (``R-H + H <=> R + H2``) is resolvable + while a hydrogen bound to something else stays excluded as a spectator. + + ``H`` is the hydrogen whose two partners are closest to it in sum. Every tie, both between two + candidate hydrogens and between two equidistant partners of one hydrogen, is resolved in favour + of the lowest index. + + Returns: + Optional[Dict[str, int]]: ``{'H': int, 'A': int, 'B': int}``, or ``None``. + """ + symbols = xyz.get('symbols') if isinstance(xyz, dict) else None + if not symbols: + return None + dmat = xyz_to_dmat(xyz) + if dmat is None: + return None + + hydrogen_indices = [i for i, symbol in enumerate(symbols) if symbol == 'H'] + min_distance = float('inf') + selected_hydrogen = None + selected_partners = None + for hydrogen_index in hydrogen_indices: + partners = sorted( + (atom for atom in range(len(symbols)) + if atom != hydrogen_index + and (symbols[atom] != 'H' + or _is_free_h_atom(symbols=symbols, + dmat=dmat, + index=atom, + transferred_h_index=hydrogen_index, + ))), + key=lambda atom: dmat[hydrogen_index][atom], + )[:2] + if len(partners) < 2: + continue + distances = dmat[hydrogen_index][partners[0]] + dmat[hydrogen_index][partners[1]] + if distances < min_distance: + min_distance = distances + selected_hydrogen = hydrogen_index + selected_partners = partners + + if selected_hydrogen is not None and selected_partners is not None: + return {'H': selected_hydrogen, 'A': selected_partners[0], 'B': selected_partners[1]} + + logger.warning('No valid hydrogen atom found for CREST H-abstraction atoms.') + return None diff --git a/arc/job/adapters/ts/xy_addition.py b/arc/job/adapters/ts/xy_addition.py new file mode 100644 index 0000000000..bb3c468c9c --- /dev/null +++ b/arc/job/adapters/ts/xy_addition.py @@ -0,0 +1,305 @@ +""" +TS-guess seed builder for the ``XY_Addition_MultipleBond`` reaction family. + +The family adds an X-Y bond across a multiple bond. Its RMG recipe is:: + + ['BREAK_BOND', '*3', 1, '*4'], # break the X-Y bond + ['CHANGE_BOND', '*1', -1, '*2'], # reduce the multiple bond + ['FORM_BOND', '*1', 1, '*3'], # form *1-X + ['FORM_BOND', '*2', 1, '*4'], # form *2-Y + +so the transition state is a 4-center arrangement: X (``*3``) approaching one end of the +multiple bond (``*1``) and Y (``*4``) approaching the other (``*2``), while the X-Y bond +(``*3``-``*4``) breaks. This module builds that 4-center geometry from the reactant +geometries and the family's atom labels, to seed a downstream TS search (e.g. CREST +refinement followed by a saddle-point optimization). + +The builder handles the bimolecular case (the multiple bond and the X-Y group are on +separate reactants); other topologies, including reactions with more or fewer than two +reactants, are skipped. + +Seed distances are set as multiples of the equilibrium single-bond length of the atom pair +involved rather than as absolute lengths, so that the seed is TS-like for any member of the +family (``*3`` is H or a halogen and ``*4`` is always a halogen, so their equilibrium bond +lengths span roughly 1.1-1.9 A). The module constants are: + +* ``FORMING_X_FACTOR`` -- multiplier for the forming ``*1``...``*3`` bond. +* ``FORMING_Y_FACTOR`` -- multiplier for the forming ``*2``...``*4`` bond. +* ``BREAKING_XY_FACTOR`` -- multiplier for the breaking ``*3``-``*4`` bond. It is larger than + one so the X-Y fragment is stretched past its equilibrium length; an unstretched fragment + yields a van der Waals complex rather than a saddle-point candidate, since bond order is + created at ``*1`` and ``*2`` without any being released at ``*3``-``*4``. +* ``LOCAL_PLANE_RADIUS`` -- radius in Angstrom around the multiple-bond midpoint within which + atoms (the bond atoms and their substituents) define the local plane of the multiple bond. + +The two forming-bond multipliers differ from each other, which relies on the family's group +definition: ``*4`` is restricted to ``[F1s,Cl1s,Br1s]`` while ``*3`` is ``[H,F1s,Cl1s,Br1s]``, +so RMG labels the hydrogen of an H-X fragment ``*3``. Their values reproduce the concerted +C2H4 + HCl four-centre transition state, in which the forming C...H bond is much further +advanced than the forming C...halogen bond. When both X and Y are halogens RMG enumerates both +``*3``/``*4`` assignments and a seed is built for each. +""" + +from typing import TYPE_CHECKING + +import numpy as np + +from arc.common import get_logger, get_single_bond_length +from arc.species.species import colliding_atoms + +if TYPE_CHECKING: + from arc.reaction import ARCReaction + +logger = get_logger() + +FORMING_X_FACTOR = 1.31 +FORMING_Y_FACTOR = 1.50 +BREAKING_XY_FACTOR = 1.42 +LOCAL_PLANE_RADIUS = 2.6 + + +def xy_addition(reaction: 'ARCReaction') -> list[dict]: + """ + Generate 4-center TS-guess seeds for an ``XY_Addition_MultipleBond`` reaction. + + One seed is built per product dictionary that carries the four family labels + (``*1``, ``*2``, ``*3``, ``*4``) and whose multiple bond and X-Y group sit on + different reactants (the bimolecular case). + + Reactions with a number of reactants other than two return no seeds: the seed geometry is + reassembled from exactly two fragments, so any additional reactant would have its symbols + and coordinates drawn from the X-Y fragment. + + Args: + reaction (ARCReaction): The reaction. Must have ``product_dicts`` populated + (each with an ``r_label_map``) and reactant geometries. + + Returns: + list[dict]: Seed entries with Cartesian coordinates and the explicit family-label + atom mapping used to build each geometry. + """ + seeds: list[dict] = list() + reactants, _ = reaction.get_reactants_and_products(return_copies=True) + if len(reactants) != 2: + logger.info(f'The XY-addition seed builder only handles bimolecular reactions, ' + f'got {len(reactants)} reactant(s) for {reaction.label}.') + return seeds + lengths = [spc.number_of_atoms for spc in reactants] + offsets = [sum(lengths[:i]) for i in range(len(reactants))] + total_atoms = sum(lengths) + + def which_reactant(global_index: int) -> tuple[int | None, int | None]: + """Map a global reactant atom index to (reactant index, local atom index).""" + for reactant_index, offset in enumerate(offsets): + if offset <= global_index < offset + lengths[reactant_index]: + return reactant_index, global_index - offset + return None, None + + for product_dict in reaction.product_dicts: + r_label_map = product_dict.get('r_label_map', dict()) + if not all(label in r_label_map for label in ('*1', '*2', '*3', '*4')): + continue + (r_mb, i1), (r_mb2, i2) = which_reactant(r_label_map['*1']), which_reactant(r_label_map['*2']) + (r_xy, i3), (r_xy2, i4) = which_reactant(r_label_map['*3']), which_reactant(r_label_map['*4']) + # The multiple bond (*1, *2) must share one reactant and the X-Y group (*3, *4) + # the other. Only the bimolecular case is handled here. + if r_mb is None or r_mb != r_mb2 or r_xy != r_xy2 or r_mb == r_xy: + continue + mb_xyz, xy_xyz = reactants[r_mb].get_xyz(), reactants[r_xy].get_xyz() + mb_coords = np.array(mb_xyz['coords'], dtype=float) + xy_coords = np.array(xy_xyz['coords'], dtype=float) + placed_xy_coords = _build_4_center_geometry(mb_coords, i1, i2, xy_coords, i3, i4, + mb_symbols=tuple(mb_xyz['symbols']), + xy_symbols=tuple(xy_xyz['symbols'])) + if placed_xy_coords is None: + continue + # Reassemble a single geometry in the reaction's reactant (global) atom order. + symbols, isotopes, coords = list(), list(), list() + for global_index in range(total_atoms): + reactant_index, local_index = which_reactant(global_index) + source_xyz = mb_xyz if reactant_index == r_mb else xy_xyz + source_coords = mb_coords if reactant_index == r_mb else placed_xy_coords + symbols.append(source_xyz['symbols'][local_index]) + isotopes.append(source_xyz['isotopes'][local_index]) + coords.append(tuple(float(v) for v in source_coords[local_index])) + seed_xyz = {'symbols': tuple(symbols), + 'isotopes': tuple(isotopes), + 'coords': tuple(coords)} + if colliding_atoms(seed_xyz): + # The X-Y fragment is placed rigidly, so it can land on a substituent. Every other seed + # builder rejects colliding geometries rather than sending them into a TS optimization. + logger.info(f'Discarding an XY-addition seed for {reaction.label}: it has colliding atoms.') + continue + seeds.append({'xyz': seed_xyz, + 'method': 'Heuristics-XY', + 'metadata': { + 'reactive_atoms': { + label: r_label_map[label] for label in ('*1', '*2', '*3', '*4') + }, + }}) + return seeds + + +def _build_4_center_geometry(mb_coords: np.ndarray, + i1: int, + i2: int, + xy_coords: np.ndarray, + i3: int, + i4: int, + mb_symbols: tuple, + xy_symbols: tuple, + ) -> np.ndarray | None: + """ + Position the X-Y fragment over the multiple bond in a 4-center arrangement. + + The multiple-bond fragment is kept fixed. The X-Y fragment is stretched onto its breaking-bond + distance, then translated and rotated so that *3 and *4 sit over *1 and *2 at their forming-bond + distances, approaching the pi face. All three distances scale with the equilibrium single-bond + length of the pair involved, so the seed is TS-like for any member of the family rather than only + for those whose atoms happen to match a fixed constant. + + ``FORMING_X_FACTOR`` is applied to ``*1``...``*3`` and ``FORMING_Y_FACTOR`` to ``*2``...``*4``, + which assumes the family's label convention that ``*4`` is a halogen while ``*3`` may be the + hydrogen of an H-X fragment. + + Args: + mb_coords (np.ndarray): Coordinates of the multiple-bond reactant. + i1, i2 (int): Local indices of the multiple-bond atoms (``*1``, ``*2``). + xy_coords (np.ndarray): Coordinates of the X-Y reactant. + i3, i4 (int): Local indices of the X-Y atoms (``*3`` = X, ``*4`` = Y). + mb_symbols (tuple): Chemical element symbols of the multiple-bond reactant. + xy_symbols (tuple): Chemical element symbols of the X-Y reactant. + + Returns: + np.ndarray | None: The transformed X-Y coordinates, or ``None`` if the multiple bond is + degenerate, if the X-Y bond is degenerate, or if the three seed + distances cannot close a ring. + """ + p1, p2 = mb_coords[i1], mb_coords[i2] + bond = p2 - p1 + bond_length = np.linalg.norm(bond) + if bond_length < 1e-3: + return None + bond_axis = bond / bond_length + midpoint = (p1 + p2) / 2 + # Approach normal: perpendicular to the local plane of the multiple bond and its substituents, + # i.e. the pi face. Projecting a vector towards the centroid perpendicular to the bond axis + # instead leaves the normal *inside* the molecular plane for any planar alkene, which sends the + # X-Y fragment edge-on into the substituents. + local = np.array([coord for coord in mb_coords + if np.linalg.norm(coord - midpoint) < LOCAL_PLANE_RADIUS]) + normal = np.zeros(3) + if len(local) >= 3: + # The smallest right singular vector of the centred local coordinates is the direction of + # least variance, i.e. the normal of the plane those atoms lie in. + normal = np.linalg.svd(local - local.mean(axis=0))[2][-1] + normal = normal - np.dot(normal, bond_axis) * bond_axis + if np.linalg.norm(normal) < 1e-3: + # A linear or two-atom fragment has no plane; any direction perpendicular to the bond works. + normal = np.cross(bond_axis, np.array([0.0, 0.0, 1.0])) + if np.linalg.norm(normal) < 1e-3: + normal = np.cross(bond_axis, np.array([0.0, 1.0, 0.0])) + normal = normal / np.linalg.norm(normal) + d_13 = FORMING_X_FACTOR * get_single_bond_length(mb_symbols[i1], xy_symbols[i3]) + d_24 = FORMING_Y_FACTOR * get_single_bond_length(mb_symbols[i2], xy_symbols[i4]) + d_34 = BREAKING_XY_FACTOR * get_single_bond_length(xy_symbols[i3], xy_symbols[i4]) + # Both faces are chemically equivalent; approach the less hindered one. + others = np.array([coord for i, coord in enumerate(mb_coords) if i not in (i1, i2)]) + if len(others): + clearances = [min(np.linalg.norm(others - (p1 + d_13 * sign * normal), axis=1)) + for sign in (1.0, -1.0)] + if clearances[1] > clearances[0]: + normal = -normal + placed = _solve_ring_positions(p1=p1, p2=p2, bond_axis=bond_axis, normal=normal, + d_13=d_13, d_24=d_24, d_34=d_34) + if placed is None: + return None + target_x, target_y = placed + # Stretch the X-Y fragment onto the breaking-bond distance before placing it, then translate and + # rotate it rigidly. Members of this family are diatomic, so the stretch is exact for them. + xy_axis = xy_coords[i4] - xy_coords[i3] + xy_length = np.linalg.norm(xy_axis) + if xy_length < 1e-3: + return None + stretched = xy_coords + (d_34 - xy_length) * (xy_axis / xy_length) * ( + np.arange(len(xy_coords)) != i3).astype(float)[:, None] + translated = stretched - stretched[i3] + target_x + rotation = _rotation_matrix_between(translated[i4] - target_x, target_y - target_x) + return (rotation @ (translated - target_x).T).T + target_x + + +def _solve_ring_positions(p1: np.ndarray, + p2: np.ndarray, + bond_axis: np.ndarray, + normal: np.ndarray, + d_13: float, + d_24: float, + d_34: float, + ) -> tuple[np.ndarray, np.ndarray] | None: + """ + Place *3 and *4 in the ring plane so all three seed distances are satisfied simultaneously. + + Works in the 2D frame spanned by the multiple-bond axis and the approach normal, with *1 at the + origin. *3 sits directly over *1; *4 is the intersection of the circle of radius ``d_24`` about + *2 with the circle of radius ``d_34`` about *3, taking the solution on the same face as *3. + + Args: + p1, p2 (np.ndarray): Coordinates of *1 and *2. + bond_axis (np.ndarray): Unit vector from *1 to *2. + normal (np.ndarray): Unit approach normal, perpendicular to ``bond_axis``. + d_13, d_24, d_34 (float): Target *1-*3, *2-*4 and *3-*4 distances. + + Returns: + tuple[np.ndarray, np.ndarray] | None: Coordinates of *3 and *4, or ``None`` if the three + distances cannot close a ring on this bond. + """ + length = float(np.linalg.norm(p2 - p1)) + x_uv = np.array([0.0, d_13]) + centre_distance = float(np.linalg.norm(x_uv - np.array([length, 0.0]))) + if centre_distance > d_34 + d_24 or centre_distance < abs(d_34 - d_24) or centre_distance < 1e-6: + return None + # Standard two-circle intersection, expressed along and perpendicular to the centre-centre line. + a = (d_34 ** 2 - d_24 ** 2 + centre_distance ** 2) / (2 * centre_distance) + h_squared = d_34 ** 2 - a ** 2 + if h_squared < 0: + return None + along = (np.array([length, 0.0]) - x_uv) / centre_distance + perpendicular = np.array([-along[1], along[0]]) + base = x_uv + a * along + candidates = [base + np.sqrt(h_squared) * perpendicular, base - np.sqrt(h_squared) * perpendicular] + y_uv = max(candidates, key=lambda candidate: candidate[1]) + to_3d = lambda uv: p1 + uv[0] * bond_axis + uv[1] * normal + return to_3d(x_uv), to_3d(y_uv) + + +def _rotation_matrix_between(vector_from: np.ndarray, vector_to: np.ndarray) -> np.ndarray: + """ + Return the rotation matrix that aligns ``vector_from`` onto ``vector_to`` (Rodrigues). + + Args: + vector_from (np.ndarray): The source vector. + vector_to (np.ndarray): The target vector. + + Returns: + np.ndarray: A 3x3 rotation matrix. + """ + a = vector_from / np.linalg.norm(vector_from) + b = vector_to / np.linalg.norm(vector_to) + axis = np.cross(a, b) + sine = np.linalg.norm(axis) + cosine = float(np.dot(a, b)) + if sine < 1e-8: + if cosine > 0: + return np.eye(3) + # Antiparallel: rotate by 180 degrees about any axis perpendicular to a. Returning -I here + # would be an improper transform (det = -1), i.e. the mirror image of the fragment. + perpendicular = np.array([1.0, 0.0, 0.0]) + if abs(float(np.dot(a, perpendicular))) > 0.9: + perpendicular = np.array([0.0, 1.0, 0.0]) + axis_180 = np.cross(a, perpendicular) + axis_180 = axis_180 / np.linalg.norm(axis_180) + return 2.0 * np.outer(axis_180, axis_180) - np.eye(3) + skew = np.array([[0.0, -axis[2], axis[1]], + [axis[2], 0.0, -axis[0]], + [-axis[1], axis[0], 0.0]]) + return np.eye(3) + skew + skew @ skew * ((1.0 - cosine) / (sine * sine)) diff --git a/arc/job/adapters/ts/xy_addition_test.py b/arc/job/adapters/ts/xy_addition_test.py new file mode 100644 index 0000000000..755eed8658 --- /dev/null +++ b/arc/job/adapters/ts/xy_addition_test.py @@ -0,0 +1,353 @@ +#!/usr/bin/env python3 +# encoding: utf-8 + +""" +Unit tests for the arc.job.adapters.ts.xy_addition module (XY_Addition_MultipleBond TS seed builder). +""" + +import math +import unittest +from unittest.mock import patch + +import numpy as np + +from arc.common import get_single_bond_length +from arc.job.adapters.ts.xy_addition import (BREAKING_XY_FACTOR, + FORMING_X_FACTOR, + FORMING_Y_FACTOR, + LOCAL_PLANE_RADIUS, + _build_4_center_geometry, + _rotation_matrix_between, + _solve_ring_positions, + xy_addition, + ) +from arc.job.adapters.ts.seed_hub import get_ts_seeds, get_wrapper_constraints +from arc.reaction import ARCReaction +from arc.species import ARCSpecies + + +ETHYLENE_XYZ = ('C 0 0 0.667\nC 0 0 -0.667\nH 0 0.921 1.232\n' + 'H 0 -0.921 1.232\nH 0 0.921 -1.232\nH 0 -0.921 -1.232') +HCL_XYZ = 'Cl 0 0 0.071\nH 0 0 -1.211' +CHLOROETHANE_XYZ = ('C 1.61 -0.36 0\nC 0.49 0.66 0\nCl -1.14 -0.15 0\nH 1.57 -0.99 -0.88\n' + 'H 2.57 0.16 0\nH 0.53 1.30 0.88\nH 0.53 1.30 -0.88\nH -0.34 -0.5 0') +HCL_COORDS = np.array([(0.0, 0.0, 0.071), (0.0, 0.0, -1.211)], dtype=float) +HCL_SYMBOLS = ('Cl', 'H') +XY_FAMILY = 'XY_Addition_MultipleBond' + + +def distance(coords, i: int, j: int) -> float: + """ + Return the distance between two atoms of a coordinates sequence. + + Args: + coords: A sequence of Cartesian coordinate triples. + i (int): The index of the first atom. + j (int): The index of the second atom. + + Returns: + float: The interatomic distance. + """ + return math.sqrt(sum((coords[i][k] - coords[j][k]) ** 2 for k in range(3))) + + +class TestXYAdditionSeed(unittest.TestCase): + """Tests for the XY_Addition_MultipleBond 4-center TS seed builder.""" + + @classmethod + def setUpClass(cls): + """Set up the shared C2H4 + HCl <=> CH3CH2Cl fixture.""" + cls.maxDiff = None + cls.label_maps = [ + {'*2': 0, '*1': 1, '*3': 7, '*4': 6}, + {'*2': 1, '*1': 0, '*3': 7, '*4': 6}, + ] + + def make_ethylene_hcl_reaction(self) -> ARCReaction: + """ + Build the C2H4 + HCl <=> CH3CH2Cl reaction with the family labels RMG assigns to it. + + The ``r_label_map`` values match the ones ARC's family module generates for this + reaction: the multiple-bond carbons are ``*1``/``*2`` (in both orderings), the hydrogen + of HCl is ``*3`` and its chlorine is ``*4``. + + Returns: + ARCReaction: The reaction, with ``product_dicts`` and ``family`` populated. + """ + ethylene = ARCSpecies(label='R1', smiles='C=C', multiplicity=1, xyz=ETHYLENE_XYZ) + hcl = ARCSpecies(label='R2', smiles='Cl', multiplicity=1, xyz=HCL_XYZ) + chloroethane = ARCSpecies(label='P1', smiles='CCCl', multiplicity=1, xyz=CHLOROETHANE_XYZ) + rxn = ARCReaction(r_species=[ethylene, hcl], p_species=[chloroethane]) + rxn.product_dicts = [{'family': XY_FAMILY, 'r_label_map': label_map} for label_map in self.label_maps] + rxn.family = XY_FAMILY + return rxn + + def test_xy_addition_seed_is_four_center(self): + """The seed for C=C + HCl -> CH3CH2Cl must be a 4-center arrangement: + X (H) and Y (Cl) both approaching the (former) double bond, with the H-Cl bond STRETCHED. + + A rigid, unstretched H-Cl makes the seed a van der Waals complex rather than a saddle point: + bond order is created at the two carbons without any being released at H-Cl, so the + imaginary mode is a hindered rotation of HCl over the pi cloud instead of the reaction + coordinate. The breaking bond must therefore be well beyond its equilibrium length.""" + rxn = self.make_ethylene_hcl_reaction() + seeds = xy_addition(reaction=rxn) + self.assertEqual(len(seeds), 2) + self.assertEqual(seeds[0]['method'], 'Heuristics-XY') + self.assertEqual( + [seed['metadata']['reactive_atoms'] for seed in seeds], + self.label_maps, + ) + hub_seeds = get_ts_seeds(reaction=rxn) + self.assertEqual([seed['metadata']['reactive_atoms'] for seed in hub_seeds], self.label_maps) + self.assertEqual( + [get_wrapper_constraints('crest', reaction=rxn, seed=seed) for seed in hub_seeds], + [ + { + 'atoms': tuple(label_map[label] for label in ('*1', '*2', '*3', '*4')), + 'distance_pairs': ( + (label_map['*1'], label_map['*3']), + (label_map['*2'], label_map['*4']), + (label_map['*3'], label_map['*4']), + ), + } + for label_map in self.label_maps + ], + ) + xyz = seeds[0]['xyz'] + self.assertEqual(len(xyz['symbols']), 8) + + coords = [tuple(c) for c in xyz['coords']] + + def dist(i, j): + return distance(coords, i, j) + + carbons = [i for i, s in enumerate(xyz['symbols']) if s == 'C'] + cl = [i for i, s in enumerate(xyz['symbols']) if s == 'Cl'][0] + hydrogens = [i for i, s in enumerate(xyz['symbols']) if s == 'H'] + transferring_h = min(hydrogens, key=lambda h: dist(h, cl)) + + min_cl_c = min(dist(cl, c) for c in carbons) + min_h_c = min(dist(transferring_h, c) for c in carbons) + # 4-center TS region: both Cl and the transferring H engage the carbons (forming bonds), + # while the H-Cl bond is still present (breaking). Contrast with an H-transfer saddle where + # Cl would be a spectator (> 2.9 A from every carbon). + # A reference concerted C2H4 + HCl four-centre TS has r(C-Cl) ~ 2.5-2.8 A: the forming bond + # is far from complete, but the Cl is not a spectator either. + self.assertLess(min_cl_c, 2.9, 'Cl should be forming a bond to a carbon (not a spectator)') + self.assertGreater(min_cl_c, 2.0, 'the C-Cl bond must not be near-formed in the seed') + self.assertLess(min_h_c, 1.9, 'the transferring H should be forming a bond to a carbon') + h_cl = dist(transferring_h, cl) + self.assertLess(h_cl, 2.2, 'the breaking H-Cl bond should still be present') + # 1.275 A is the equilibrium H-Cl length; an unstretched bond means a reactant complex. + self.assertGreater(h_cl, 1.40, 'the breaking H-Cl bond must be stretched well past equilibrium') + + def test_seed_orientation_follows_the_label_map(self): + """*3 must approach *1 and *4 must approach *2, at the scaled seed distances. + + Resolving the atoms through the label map (rather than by element) is what makes this + directional: a seed built with *3 and *4 - or *1 and *2 - interchanged contradicts the + CREST constraint spec that ``get_wrapper_constraints`` derives from the same label map. + """ + rxn = self.make_ethylene_hcl_reaction() + seeds = xy_addition(reaction=rxn) + self.assertEqual(len(seeds), 2) + expected_d_13 = 1.31 * get_single_bond_length('C', 'H') + expected_d_24 = 1.50 * get_single_bond_length('C', 'Cl') + expected_d_34 = 1.42 * get_single_bond_length('H', 'Cl') + for seed, label_map in zip(seeds, self.label_maps): + with self.subTest(label_map=label_map): + coords = [tuple(c) for c in seed['xyz']['coords']] + i1, i2 = label_map['*1'], label_map['*2'] + i3, i4 = label_map['*3'], label_map['*4'] + self.assertEqual(seed['xyz']['symbols'][i3], 'H') + self.assertEqual(seed['xyz']['symbols'][i4], 'Cl') + self.assertLess(distance(coords, i1, i3), distance(coords, i2, i3), + 'X (*3) must approach *1, not *2') + self.assertLess(distance(coords, i2, i4), distance(coords, i1, i4), + 'Y (*4) must approach *2, not *1') + self.assertAlmostEqual(distance(coords, i1, i3), expected_d_13, delta=0.05) + self.assertAlmostEqual(distance(coords, i2, i4), expected_d_24, delta=0.05) + self.assertAlmostEqual(distance(coords, i3, i4), expected_d_34, delta=0.05) + + def test_seed_scaling_factors(self): + """The seed scaling factors are pinned; the geometry tests are written against these values.""" + self.assertAlmostEqual(FORMING_X_FACTOR, 1.31, delta=1e-9) + self.assertAlmostEqual(FORMING_Y_FACTOR, 1.50, delta=1e-9) + self.assertAlmostEqual(BREAKING_XY_FACTOR, 1.42, delta=1e-9) + self.assertAlmostEqual(LOCAL_PLANE_RADIUS, 2.6, delta=1e-9) + + def test_approach_is_perpendicular_to_the_local_plane(self): + """X and Y must approach the pi face, i.e. along the normal of the local plane. + + The fixture is an ethylene rotated by 45 degrees about its C=C axis, so the plane normal + is not aligned with any Cartesian axis and a fallback direction cannot pass by accident. + """ + cos45 = math.sqrt(0.5) + mb_coords = np.array([(0.0, 0.0, 0.667), + (0.0, 0.0, -0.667), + (-0.921 * cos45, 0.921 * cos45, 1.232), + (0.921 * cos45, -0.921 * cos45, 1.232), + (-0.921 * cos45, 0.921 * cos45, -1.232), + (0.921 * cos45, -0.921 * cos45, -1.232)], dtype=float) + placed = _build_4_center_geometry(mb_coords, 0, 1, HCL_COORDS, 1, 0, + mb_symbols=('C', 'C', 'H', 'H', 'H', 'H'), + xy_symbols=HCL_SYMBOLS) + self.assertIsNotNone(placed) + approach = placed[1] - mb_coords[0] + approach = approach / np.linalg.norm(approach) + plane_normal = np.array([cos45, cos45, 0.0]) + self.assertGreater(abs(float(np.dot(approach, plane_normal))), 0.95, + 'the approach direction must be the local plane normal') + + def test_approach_is_orthogonal_to_the_bond_axis(self): + """The *1...*3 approach vector must be perpendicular to the multiple-bond axis. + + The fixture is a synthetic non-planar fragment whose least-variance direction coincides + with the multiple-bond axis, so a normal that is not orthogonalized against that axis + would place X along the bond instead of over it. + """ + mb_coords = np.array([(0.0, 0.0, 0.667), + (0.0, 0.0, -0.667), + (1.5, 0.0, 0.0), + (-1.5, 0.0, 0.0), + (0.0, 1.5, 0.0), + (0.0, -1.5, 0.0)], dtype=float) + placed = _build_4_center_geometry(mb_coords, 0, 1, HCL_COORDS, 1, 0, + mb_symbols=('C', 'C', 'H', 'H', 'H', 'H'), + xy_symbols=HCL_SYMBOLS) + self.assertIsNotNone(placed) + bond_axis = mb_coords[1] - mb_coords[0] + bond_axis = bond_axis / np.linalg.norm(bond_axis) + self.assertAlmostEqual(float(np.dot(placed[1] - mb_coords[0], bond_axis)), 0.0, delta=0.05) + + def test_the_less_hindered_face_is_selected(self): + """The X-Y fragment must approach the face that is not blocked by a substituent. + + The same fragment is placed twice with the blocking atom moved from one face to the + other; X must land on the opposite side each time. + """ + placed_x = list() + for blocker_sign in (1.0, -1.0): + mb_coords = np.array([(0.0, 0.0, 0.667), + (0.0, 0.0, -0.667), + (0.0, 1.6, 1.6), + (0.0, -1.6, -1.6), + (blocker_sign * 1.4, 0.0, 2.6)], dtype=float) + placed = _build_4_center_geometry(mb_coords, 0, 1, HCL_COORDS, 1, 0, + mb_symbols=('C', 'C', 'H', 'H', 'H'), + xy_symbols=HCL_SYMBOLS) + self.assertIsNotNone(placed) + placed_x.append(float(placed[1][0])) + self.assertLess(placed_x[0], 0.0, 'X must avoid a blocker on the +x face') + self.assertGreater(placed_x[1], 0.0, 'X must avoid a blocker on the -x face') + + def test_rotation_matrix_between_is_a_proper_rotation(self): + """``_rotation_matrix_between`` must never return an improper (mirroring) transform. + + The antiparallel case is the one at risk: -I maps the vector correctly but has a + determinant of -1, which inverts the chirality of the fragment being placed. + """ + antiparallel = _rotation_matrix_between(np.array([0.0, 0.0, 1.0]), np.array([0.0, 0.0, -1.0])) + self.assertAlmostEqual(float(np.linalg.det(antiparallel)), 1.0, delta=1e-8) + np.testing.assert_allclose(antiparallel @ np.array([0.0, 0.0, 1.0]), + np.array([0.0, 0.0, -1.0]), atol=1e-8) + np.testing.assert_allclose(antiparallel.T @ antiparallel, np.eye(3), atol=1e-8) + antiparallel_x = _rotation_matrix_between(np.array([1.0, 0.0, 0.0]), np.array([-2.0, 0.0, 0.0])) + self.assertAlmostEqual(float(np.linalg.det(antiparallel_x)), 1.0, delta=1e-8) + np.testing.assert_allclose(antiparallel_x @ np.array([1.0, 0.0, 0.0]), + np.array([-1.0, 0.0, 0.0]), atol=1e-8) + parallel = _rotation_matrix_between(np.array([0.0, 0.0, 1.0]), np.array([0.0, 0.0, 3.0])) + np.testing.assert_allclose(parallel, np.eye(3), atol=1e-8) + general = _rotation_matrix_between(np.array([1.0, 0.0, 0.0]), np.array([0.0, 2.0, 0.0])) + self.assertAlmostEqual(float(np.linalg.det(general)), 1.0, delta=1e-8) + np.testing.assert_allclose(general @ np.array([1.0, 0.0, 0.0]), + np.array([0.0, 1.0, 0.0]), atol=1e-8) + + def test_no_seeds_when_all_labels_are_on_one_reactant(self): + """A unimolecular arrangement of the four labels yields no seed.""" + rxn = self.make_ethylene_hcl_reaction() + rxn.product_dicts = [{'family': XY_FAMILY, 'r_label_map': {'*1': 0, '*2': 1, '*3': 2, '*4': 3}}] + self.assertEqual(xy_addition(reaction=rxn), list()) + + def test_no_seeds_when_labels_are_incomplete(self): + """A product dict missing one of the four family labels yields no seed.""" + rxn = self.make_ethylene_hcl_reaction() + rxn.product_dicts = [{'family': XY_FAMILY, 'r_label_map': {'*1': 1, '*2': 0, '*3': 7}}] + self.assertEqual(xy_addition(reaction=rxn), list()) + + def test_no_seeds_for_more_than_two_reactants(self): + """Only bimolecular reactions are assembled. + + The reassembly loop resolves every atom that is not on the multiple-bond reactant from + the X-Y fragment, so a third reactant would take its symbols and coordinates from the + wrong fragment. + """ + ethylene = ARCSpecies(label='R1', smiles='C=C', multiplicity=1, xyz=ETHYLENE_XYZ) + hcl = ARCSpecies(label='R2', smiles='Cl', multiplicity=1, xyz=HCL_XYZ) + water = ARCSpecies(label='R3', smiles='O', multiplicity=1, + xyz='O 0 0 0.119\nH 0 0.763 -0.477\nH 0 -0.763 -0.477') + chloroethane = ARCSpecies(label='P1', smiles='CCCl', multiplicity=1, xyz=CHLOROETHANE_XYZ) + water_p = ARCSpecies(label='P2', smiles='O', multiplicity=1, + xyz='O 0 0 0.119\nH 0 0.763 -0.477\nH 0 -0.763 -0.477') + rxn = ARCReaction(r_species=[ethylene, hcl, water], p_species=[chloroethane, water_p]) + rxn.product_dicts = [{'family': XY_FAMILY, 'r_label_map': label_map} for label_map in self.label_maps] + rxn.family = XY_FAMILY + self.assertEqual(xy_addition(reaction=rxn), list()) + + def test_colliding_seeds_are_discarded(self): + """A rigidly placed X-Y fragment that lands on another atom is not returned.""" + rxn = self.make_ethylene_hcl_reaction() + with patch('arc.job.adapters.ts.xy_addition.colliding_atoms', return_value=True): + self.assertEqual(xy_addition(reaction=rxn), list()) + self.assertEqual(len(xy_addition(reaction=rxn)), 2) + + def test_solve_ring_positions_returns_none_when_infeasible(self): + """Distances that cannot close the four-membered ring give ``None``, not NaN coordinates.""" + p1, p2 = np.zeros(3), np.array([1.33, 0.0, 0.0]) + bond_axis, normal = np.array([1.0, 0.0, 0.0]), np.array([0.0, 1.0, 0.0]) + for d_13, d_24, d_34 in ((1.43, 2.65, 0.2), (1.43, 0.2, 0.2), (1.43, 12.0, 1.8)): + with self.subTest(d_13=d_13, d_24=d_24, d_34=d_34): + self.assertIsNone(_solve_ring_positions(p1=p1, p2=p2, bond_axis=bond_axis, normal=normal, + d_13=d_13, d_24=d_24, d_34=d_34)) + feasible = _solve_ring_positions(p1=p1, p2=p2, bond_axis=bond_axis, normal=normal, + d_13=1.43, d_24=2.65, d_34=1.8) + self.assertIsNotNone(feasible) + for position in feasible: + self.assertTrue(np.all(np.isfinite(position))) + + def test_build_4_center_geometry_returns_none_for_degenerate_bonds(self): + """Coincident multiple-bond atoms or coincident X-Y atoms give ``None``.""" + degenerate_mb = np.array([(0.0, 0.0, 0.0), + (0.0, 0.0, 0.0), + (0.0, 0.921, 1.232), + (0.0, -0.921, 1.232)], dtype=float) + self.assertIsNone(_build_4_center_geometry(degenerate_mb, 0, 1, HCL_COORDS, 1, 0, + mb_symbols=('C', 'C', 'H', 'H'), + xy_symbols=HCL_SYMBOLS)) + mb_coords = np.array([(0.0, 0.0, 0.667), + (0.0, 0.0, -0.667), + (0.0, 0.921, 1.232), + (0.0, -0.921, 1.232), + (0.0, 0.921, -1.232), + (0.0, -0.921, -1.232)], dtype=float) + degenerate_xy = np.array([(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)], dtype=float) + self.assertIsNone(_build_4_center_geometry(mb_coords, 0, 1, degenerate_xy, 1, 0, + mb_symbols=('C', 'C', 'H', 'H', 'H', 'H'), + xy_symbols=HCL_SYMBOLS)) + + def test_build_4_center_geometry_returns_none_when_the_ring_cannot_close(self): + """An X-Y fragment that cannot reach both multiple-bond atoms gives ``None``.""" + mb_coords = np.array([(0.0, 0.0, 6.0), + (0.0, 0.0, -6.0), + (0.0, 0.921, 7.0), + (0.0, -0.921, 7.0), + (0.0, 0.921, -7.0), + (0.0, -0.921, -7.0)], dtype=float) + self.assertIsNone(_build_4_center_geometry(mb_coords, 0, 1, HCL_COORDS, 1, 0, + mb_symbols=('C', 'C', 'H', 'H', 'H', 'H'), + xy_symbols=HCL_SYMBOLS)) + + +if __name__ == '__main__': + unittest.main() diff --git a/arc/scheduler.py b/arc/scheduler.py index f606976396..a69244839a 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -2427,6 +2427,11 @@ def determine_most_likely_ts_conformer(self, label: str): Determine the most likely TS conformer. Save the resulting xyz as the ``.initial_xyz`` attribute of the TS Species. + A successful guess credits its method and every method listed in its ``method_sources`` + to the species' ``successful_methods``, deduplicated case-insensitively. A guess's method + is appended to ``unsuccessful_methods`` only if it was not credited, compared + case-insensitively as well. + Args: label (str): The TS species label. """ @@ -2437,9 +2442,14 @@ def determine_most_likely_ts_conformer(self, label: str): # Only run this block once, not every time a TS is selecting a different guess. for tsg in self.species_dict[label].ts_guesses: if tsg.success: - self.species_dict[label].successful_methods.append(tsg.method) + recorded = [method.lower() for method in self.species_dict[label].successful_methods] + for method in [tsg.method] + list(tsg.method_sources or []): + if method.lower() not in recorded: + self.species_dict[label].successful_methods.append(method) + recorded.append(method.lower()) + successful_lower = [method.lower() for method in self.species_dict[label].successful_methods] for tsg in self.species_dict[label].ts_guesses: - if tsg.method not in self.species_dict[label].successful_methods: + if tsg.method.lower() not in successful_lower: self.species_dict[label].unsuccessful_methods.append(tsg.method) message = f'\nAll TS guesses for {label} terminated.' if self.species_dict[label].successful_methods and not self.species_dict[label].unsuccessful_methods: diff --git a/arc/scheduler_test.py b/arc/scheduler_test.py index 7803e010c2..4cf16a3ba8 100644 --- a/arc/scheduler_test.py +++ b/arc/scheduler_test.py @@ -1112,6 +1112,106 @@ def test_run_ts_conformer_jobs_single_success_provenance(self, mock_run_opt): self.assertIsNone(failed.energy) self.assertEqual(sched.output['TS_single']['paths']['neb'], good.log_path) + def _make_ts_scheduler(self, label: str, project: str, folder: str, ts_guesses: list) -> Scheduler: + """ + Create a Scheduler holding a single TS species with the given TS guesses. + + Args: + label (str): The TS species label. + project (str): The project name. + folder (str): The name of the project folder under ARC's Projects directory. + ts_guesses (list): The TSGuess object instances to assign to the TS species. + + Returns: + Scheduler: The scheduler instance. + """ + ts_spc = ARCSpecies(label=label, is_ts=True, multiplicity=1, charge=0, compute_thermo=False) + ts_spc.ts_guesses = ts_guesses + project_directory = os.path.join(ARC_PATH, 'Projects', folder) + self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True) + return Scheduler(project=project, ess_settings=self.ess_settings, + species_list=[ts_spc], + opt_level=Level(repr=default_levels_of_theory['opt']), + freq_level=Level(repr=default_levels_of_theory['freq']), + sp_level=Level(repr=default_levels_of_theory['sp']), + ts_guess_level=Level(repr=default_levels_of_theory['ts_guesses']), + project_directory=project_directory, + testing=True, + job_types=self.job_types1, + ) + + def test_determine_most_likely_ts_conformer_credits_all_method_sources(self): + """Test that every method source of a successful TS guess is reported as a successful method.""" + xyz_1 = str_to_xyz("""N 0.91779059 0.51946178 0.00000000 + H 1.81402049 1.03819414 0.00000000 + H 0.00000000 0.00000000 0.00000000 + H 0.91779059 1.22790192 0.72426890""") + xyz_2 = str_to_xyz("""N 0.92779059 0.52946178 0.00000000 + H 1.84402049 1.05819414 0.00000000 + H 0.00000000 0.00000000 0.00000000 + H 0.93779059 1.25790192 0.75426890""") + xyz_3 = str_to_xyz("""N 0.95779059 0.55946178 0.00000000 + H 1.88402049 1.09819414 0.00000000 + H 0.00000000 0.00000000 0.00000000 + H 0.96779059 1.28790192 0.78426890""") + merged = TSGuess(index=0, method='heuristics', success=True, energy=100.0, xyz=xyz_1, + execution_time='0:00:01') + merged.method_sources = ['heuristics', 'crest'] + merged.opt_xyz = xyz_1 + merged.imaginary_freqs = [-500.0] + sourceless = TSGuess(index=1, method='xtb_gsm', success=True, energy=120.0, xyz=xyz_2, + execution_time='0:00:01') + sourceless.method_sources = None + sourceless.opt_xyz = xyz_2 + sourceless.imaginary_freqs = [-400.0] + failed = TSGuess(index=2, method='AutoTST', success=False, xyz=xyz_3, execution_time='0:00:01') + + label = 'TS_method_sources' + sched = self._make_ts_scheduler(label=label, project='test_ts_method_sources', + folder='arc_project_for_testing_delete_after_usage_method_sources', + ts_guesses=[merged, sourceless, failed]) + sched.determine_most_likely_ts_conformer(label=label) + + successful_methods = sched.species_dict[label].successful_methods + unsuccessful_methods = sched.species_dict[label].unsuccessful_methods + self.assertEqual(successful_methods, ['heuristics', 'crest', 'xtb_gsm']) + self.assertEqual(unsuccessful_methods, ['autotst']) + self.assertNotIn('crest', unsuccessful_methods) + self.assertNotIn('heuristics', unsuccessful_methods) + self.assertNotIn('xtb_gsm', unsuccessful_methods) + + def test_determine_most_likely_ts_conformer_dedups_methods_case_insensitively(self): + """Test that a method credited under different letter cases is reported once and never as unsuccessful.""" + xyz_1 = str_to_xyz("""N 0.91779059 0.51946178 0.00000000 + H 1.81402049 1.03819414 0.00000000 + H 0.00000000 0.00000000 0.00000000 + H 0.91779059 1.22790192 0.72426890""") + xyz_2 = str_to_xyz("""N 0.92779059 0.52946178 0.00000000 + H 1.84402049 1.05819414 0.00000000 + H 0.00000000 0.00000000 0.00000000 + H 0.93779059 1.25790192 0.75426890""") + upper = TSGuess(index=0, method='crest', success=True, energy=100.0, xyz=xyz_1, + execution_time='0:00:01') + upper.method = 'CREST' + upper.opt_xyz = xyz_1 + upper.imaginary_freqs = [-500.0] + lower = TSGuess(index=1, method='crest', success=True, energy=120.0, xyz=xyz_2, + execution_time='0:00:01') + lower.opt_xyz = xyz_2 + lower.imaginary_freqs = [-400.0] + + label = 'TS_case_dedup' + sched = self._make_ts_scheduler(label=label, project='test_ts_case_dedup', + folder='arc_project_for_testing_delete_after_usage_case_dedup', + ts_guesses=[upper, lower]) + sched.determine_most_likely_ts_conformer(label=label) + + successful_methods = sched.species_dict[label].successful_methods + unsuccessful_methods = sched.species_dict[label].unsuccessful_methods + self.assertEqual(len(successful_methods), 1) + self.assertEqual(successful_methods[0].lower(), 'crest') + self.assertEqual(unsuccessful_methods, list()) + @patch('arc.scheduler.Scheduler.run_opt_job') def test_switch_ts_cleanup(self, mock_run_opt): """Test that switch_ts resets job_types, convergence, cleans up IRC species, and clears pending pipes.""" diff --git a/arc/settings/crest.py b/arc/settings/crest.py new file mode 100644 index 0000000000..9a5abad87e --- /dev/null +++ b/arc/settings/crest.py @@ -0,0 +1,190 @@ +""" +Utilities for locating CREST executables and activation commands. +""" + +import functools +import os +import re +import shutil +import sys +from typing import Iterator, Optional, Tuple + + +CREST_ENV_NAME = "crest_env" +STANDALONE_DIR_ENV_VAR = "ARC_CREST_STANDALONE_DIR" +DEFAULT_STANDALONE_DIR = "/Local/ce_dana" + + +def parse_version(folder_name: str) -> Tuple[int, int, int]: + """ + Parse a version from a folder name. + + Supports patterns such as ``3.0.2``, ``v212``, ``2.1``, ``2``. + A three-digit run such as ``v212`` is read as ``(2, 1, 2)``. + + Args: + folder_name (str): The folder name to parse. + + Returns: + Tuple[int, int, int]: The major, minor and patch version numbers, ``(0, 0, 0)`` if no version was found. + """ + version_regex = re.compile(r"(?:v?(\d+)(?:\.(\d+))?(?:\.(\d+))?)", re.IGNORECASE) + match = version_regex.search(folder_name) + if not match: + return 0, 0, 0 + + major = int(match.group(1)) if match.group(1) else 0 + minor = int(match.group(2)) if match.group(2) else 0 + patch = int(match.group(3)) if match.group(3) else 0 + + if major >= 100 and match.group(2) is None and match.group(3) is None: + s = str(major).rjust(3, "0") + major, minor, patch = int(s[0]), int(s[1]), int(s[2]) + + return major, minor, patch + + +def find_highest_version_in_directory(directory: str, name_contains: str) -> Optional[str]: + """ + Find the ``crest`` executable under the highest-version matching subdirectory. + + Args: + directory (str): The directory to search in. + name_contains (str): A substring which a subdirectory name must contain (case-insensitively). + + Returns: + Optional[str]: The path to the executable, ``None`` if the directory is missing, + unreadable, or holds no matching executable. + """ + if not directory or not os.path.exists(directory): + return None + + highest_version_path = None + highest_version = () + try: + folders = os.listdir(directory) + except OSError: + return None + for folder in folders: + file_path = os.path.join(directory, folder) + if name_contains.lower() in folder.lower() and os.path.isdir(file_path): + crest_path = os.path.join(file_path, "crest") + if os.path.isfile(crest_path) and os.access(crest_path, os.X_OK): + version = parse_version(folder) + if highest_version == () or version > highest_version: + highest_version = version + highest_version_path = crest_path + return highest_version_path + + +def _iter_ancestor_dirs(directory: str) -> Iterator[str]: + """ + Iterate over a directory and each of its ancestors up to the filesystem root. + + Args: + directory (str): The directory to start from. + + Yields: + str: The next directory, starting with ``directory`` itself. + """ + current = directory + while current: + yield current + parent = os.path.dirname(current) + if parent == current: + break + current = parent + + +def find_env_activation_command(crest_path: str, env_name: str = CREST_ENV_NAME) -> str: + """ + Determine a shell command which activates the environment holding a crest executable. + + A command is only generated when the executable lives under an ``envs//`` directory + and an ``etc/profile.d`` activation script exists in that environment's directory or in one of + its ancestors. Otherwise an empty string is returned, and the executable is meant to be invoked + by its absolute path. Note that a user data directory such as ``~/.conda`` holds environments + but no activation script, and therefore yields an empty string. + + Args: + crest_path (str): The path to the crest executable. + env_name (str, optional): The name of the environment to activate. + + Returns: + str: A shell snippet activating the environment, an empty string if no activation script was found. + """ + env_marker = os.path.join("envs", env_name) + os.path.sep + if env_marker not in crest_path: + return "" + env_root = crest_path.split(env_marker)[0].rstrip(os.path.sep) + if not env_root: + return "" + for directory in _iter_ancestor_dirs(env_root): + micromamba_sh = os.path.join(directory, "etc", "profile.d", "micromamba.sh") + if os.path.isfile(micromamba_sh): + return f"source {micromamba_sh} && micromamba activate {env_name}" + conda_sh = os.path.join(directory, "etc", "profile.d", "conda.sh") + if os.path.isfile(conda_sh): + return f"source {conda_sh} && conda activate {env_name}" + return "" + + +def find_crest_executable() -> Tuple[Optional[str], Optional[str]]: + """ + Locate a crest executable along with the command needed to activate its environment. + + Locations are searched in this order: a standalone build under the directory named by the + ``ARC_CREST_STANDALONE_DIR`` environment variable, defaulting to ``/Local/ce_dana`` and + disabled by setting that variable to an empty string; the active Python environment and the + common conda, mamba and micromamba environment locations under the user's home directory; + and finally ``PATH``. + + Returns: + Tuple[Optional[str], Optional[str]]: The path to the crest executable and a shell snippet + activating its environment (an empty string when no activation is required), + or ``(None, None)`` when no executable was found. + """ + standalone_dir = os.getenv(STANDALONE_DIR_ENV_VAR, DEFAULT_STANDALONE_DIR) + crest_path = find_highest_version_in_directory(standalone_dir, "crest") + if crest_path and os.path.isfile(crest_path) and os.access(crest_path, os.X_OK): + return crest_path, "" + + home = os.path.expanduser("~") + potential_env_paths = [os.path.join(os.path.dirname(sys.executable), "crest")] + potential_env_paths += [os.path.join(home, root, "envs", CREST_ENV_NAME, "bin", "crest") + for root in ("anaconda3", "miniconda3", "miniforge3", ".conda", + "mambaforge", "micromamba")] + + for crest_path in potential_env_paths: + if os.path.isfile(crest_path) and os.access(crest_path, os.X_OK): + return crest_path, find_env_activation_command(crest_path) + + crest_in_path = shutil.which("crest") + if crest_in_path: + return crest_in_path, "" + + return None, None + + +@functools.lru_cache(maxsize=1) +def get_crest_paths() -> Tuple[Optional[str], Optional[str]]: + """ + Return the cached result of :func:`find_crest_executable`. + + The filesystem is only searched on the first call. Call ``get_crest_paths.cache_clear()`` + to force a new search. + + Returns: + Tuple[Optional[str], Optional[str]]: The crest executable path and its environment activation command. + """ + return find_crest_executable() + + +__all__ = [ + "DEFAULT_STANDALONE_DIR", + "parse_version", + "find_highest_version_in_directory", + "find_env_activation_command", + "find_crest_executable", + "get_crest_paths", +] diff --git a/arc/settings/crest_test.py b/arc/settings/crest_test.py new file mode 100644 index 0000000000..ac244d99d7 --- /dev/null +++ b/arc/settings/crest_test.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +# encoding: utf-8 + +""" +Unit tests for arc.settings.crest +""" + +import os +import stat +import tempfile +import unittest +from unittest.mock import patch + +from arc.settings.crest import ( + DEFAULT_STANDALONE_DIR, + find_crest_executable, + find_env_activation_command, + find_highest_version_in_directory, + get_crest_paths, + parse_version, +) + + +class TestCrestSettingsUtils(unittest.TestCase): + + def _make_executable(self, path: str): + """Create an executable stub file at ``path``.""" + os.makedirs(os.path.dirname(path), exist_ok=True) + with open(path, "w") as f: + f.write("#!/bin/bash\n") + st = os.stat(path) + os.chmod(path, st.st_mode | stat.S_IXUSR) + + def _make_activation_script(self, root: str, name: str): + """Create an ``etc/profile.d/`` activation script under ``root``.""" + script_path = os.path.join(root, "etc", "profile.d", name) + os.makedirs(os.path.dirname(script_path), exist_ok=True) + with open(script_path, "w") as f: + f.write("# activation stub\n") + return script_path + + def _isolated_env(self, standalone_dir: str = ""): + """Return a patcher pinning ARC_CREST_STANDALONE_DIR so the host filesystem cannot leak in.""" + return patch.dict(os.environ, {"ARC_CREST_STANDALONE_DIR": standalone_dir}, clear=False) + + def test_parse_version(self): + """Test parsing versions out of folder names""" + self.assertEqual(parse_version("crest-3.0.2"), (3, 0, 2)) + self.assertEqual(parse_version("v212"), (2, 1, 2)) + self.assertEqual(parse_version("version-2.1"), (2, 1, 0)) + self.assertEqual(parse_version("foo"), (0, 0, 0)) + + def test_find_highest_version_in_directory(self): + """Test that the highest-version crest build in a directory is selected""" + with tempfile.TemporaryDirectory() as td: + low = os.path.join(td, "crest-2.1") + high = os.path.join(td, "crest-3.0.2") + os.makedirs(low) + os.makedirs(high) + self._make_executable(os.path.join(low, "crest")) + self._make_executable(os.path.join(high, "crest")) + + found = find_highest_version_in_directory(td, "crest") + self.assertEqual(found, os.path.join(high, "crest")) + + def test_find_highest_version_in_directory_missing_dir(self): + """Test that a missing or unnamed directory yields None rather than raising""" + self.assertIsNone(find_highest_version_in_directory("", "crest")) + with tempfile.TemporaryDirectory() as td: + self.assertIsNone(find_highest_version_in_directory(os.path.join(td, "nope"), "crest")) + + def test_find_highest_version_in_directory_unreadable_dir(self): + """Test that an OSError from os.listdir yields None, keeping arc.settings.settings importable""" + with tempfile.TemporaryDirectory() as td: + with patch("arc.settings.crest.os.listdir", side_effect=OSError("stale file handle")): + self.assertIsNone(find_highest_version_in_directory(td, "crest")) + + def test_find_crest_executable_prefers_standalone(self): + """Test that a standalone build wins over any environment installation""" + with tempfile.TemporaryDirectory() as td: + standalone = os.path.join(td, "crest-3.0.2") + os.makedirs(standalone) + standalone_crest = os.path.join(standalone, "crest") + self._make_executable(standalone_crest) + + with self._isolated_env(standalone_dir=td): + path, env_cmd = find_crest_executable() + self.assertEqual(path, standalone_crest) + self.assertEqual(env_cmd, "") + + def test_find_crest_executable_standalone_dir_default_and_override(self): + """Test that the standalone directory defaults to /Local/ce_dana and is disabled by an empty override""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + os.makedirs(fake_home) + for standalone_dir, expected_probe in [(None, DEFAULT_STANDALONE_DIR), ("", "")]: + with patch.dict(os.environ, {}, clear=False): + if standalone_dir is None: + os.environ.pop("ARC_CREST_STANDALONE_DIR", None) + else: + os.environ["ARC_CREST_STANDALONE_DIR"] = standalone_dir + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + with patch("arc.settings.crest.find_highest_version_in_directory") as mock_find: + path, env_cmd = find_crest_executable() + mock_find.assert_called_once_with(expected_probe, "crest") + self.assertIsNone(path) + self.assertIsNone(env_cmd) + + def test_find_crest_executable_env_detection(self): + """Test detecting a conda env installation and its activation command""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + conda_root = os.path.join(fake_home, "miniforge3") + crest_path = os.path.join(conda_root, "envs", "crest_env", "bin", "crest") + self._make_executable(crest_path) + conda_sh = self._make_activation_script(conda_root, "conda.sh") + + with self._isolated_env(): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + path, env_cmd = find_crest_executable() + self.assertEqual(path, crest_path) + self.assertEqual(env_cmd, f"source {conda_sh} && conda activate crest_env") + + def test_find_crest_executable_micromamba_env_detection(self): + """Test detecting a micromamba env installation and its activation command""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + mamba_root = os.path.join(fake_home, "micromamba") + crest_path = os.path.join(mamba_root, "envs", "crest_env", "bin", "crest") + self._make_executable(crest_path) + micromamba_sh = self._make_activation_script(mamba_root, "micromamba.sh") + + with self._isolated_env(): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + path, env_cmd = find_crest_executable() + self.assertEqual(path, crest_path) + self.assertEqual(env_cmd, f"source {micromamba_sh} && micromamba activate crest_env") + + def test_find_crest_executable_dot_conda_is_not_a_conda_root(self): + """Test that ~/.conda, which holds envs but no activation script, does not yield an activation command""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + crest_path = os.path.join(fake_home, ".conda", "envs", "crest_env", "bin", "crest") + self._make_executable(crest_path) + self._make_activation_script(os.path.join(fake_home, "miniconda3"), "conda.sh") + + with self._isolated_env(): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + path, env_cmd = find_crest_executable() + self.assertEqual(path, crest_path) + self.assertEqual(env_cmd, "") + self.assertNotIn(".conda", env_cmd) + + def test_find_crest_executable_in_active_env_bin(self): + """Test that a crest in the active env bin is used as-is, without sourcing its path as a conda root""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + os.makedirs(fake_home) + active_env = os.path.join(td, "miniconda3", "envs", "arc_env") + crest_path = os.path.join(active_env, "bin", "crest") + self._make_executable(crest_path) + self._make_activation_script(os.path.join(td, "miniconda3"), "conda.sh") + + with self._isolated_env(): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(active_env, "bin", "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + path, env_cmd = find_crest_executable() + self.assertEqual(path, crest_path) + self.assertEqual(env_cmd, "") + self.assertNotIn(active_env, env_cmd) + self.assertNotIn("conda activate", env_cmd) + + def test_find_crest_executable_path_fallback(self): + """Test falling back to a crest found on PATH""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + os.makedirs(fake_home) + path_crest = os.path.join(td, "usr", "bin", "crest") + self._make_executable(path_crest) + + with self._isolated_env(): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=path_crest): + path, env_cmd = find_crest_executable() + self.assertEqual(path, path_crest) + self.assertEqual(env_cmd, "") + + def test_find_crest_executable_not_found(self): + """Test that (None, None) is returned when no crest executable exists anywhere""" + with tempfile.TemporaryDirectory() as td: + fake_home = os.path.join(td, "home") + os.makedirs(fake_home) + + with self._isolated_env(standalone_dir=os.path.join(td, "no_such_dir")): + with patch("arc.settings.crest.os.path.expanduser", return_value=fake_home): + with patch("arc.settings.crest.sys.executable", os.path.join(td, "python")): + with patch("arc.settings.crest.shutil.which", return_value=None): + path, env_cmd = find_crest_executable() + self.assertIsNone(path) + self.assertIsNone(env_cmd) + + def test_find_env_activation_command_without_env_marker(self): + """Test that an executable outside an envs/crest_env directory yields no activation command""" + self.assertEqual(find_env_activation_command("/opt/crest-3.0.2/crest"), "") + self.assertEqual(find_env_activation_command("envs/crest_env/bin/crest"), "") + + def test_get_crest_paths_is_cached(self): + """Test that the crest lookup touches the filesystem only once""" + get_crest_paths.cache_clear() + self.addCleanup(get_crest_paths.cache_clear) + with patch("arc.settings.crest.find_crest_executable", + return_value=("/opt/crest", "")) as mock_find: + first = get_crest_paths() + second = get_crest_paths() + self.assertEqual(first, ("/opt/crest", "")) + self.assertEqual(second, first) + self.assertEqual(mock_find.call_count, 1) + + +if __name__ == "__main__": + unittest.main() diff --git a/arc/settings/settings.py b/arc/settings/settings.py index 8bf3aa449e..f02d1b5386 100644 --- a/arc/settings/settings.py +++ b/arc/settings/settings.py @@ -9,6 +9,7 @@ import os import string import sys +from arc.settings.crest import get_crest_paths from arc.settings.external_paths import ( find_goflow_ckpt, @@ -116,9 +117,14 @@ supported_ess = ['cfour', 'gaussian', 'mockter', 'molpro', 'orca', 'qchem', 'terachem', 'onedmin', 'xtb', 'torchani', 'openbabel', 'ase', 'pyscf'] # TS methods to try when appropriate for a reaction (other than user guesses which are always allowed): -# Note: 'goflow' and 'rits' are intentionally NOT in the default — their envs -# (goflow_env / rits_env + pretrained checkpoints) are heavyweight, so users -# opt in explicitly via ``ts_adapters: ['goflow', ...]`` in their input.yml. +# An adapter runs only if it is both registered for the reaction's family in +# ts_adapters_by_rmg_family (arc/job/adapters/common.py) and listed here. +# 'goflow', 'rits' and 'crest' are registered per family but are not in this default list, +# so they are opt-in via ``ts_adapters: ['crest', ...]`` in the input file. +# 'heuristics' is in this default list, and it is registered for H_Abstraction, +# XY_Addition_MultipleBond, carbonyl_based_hydrolysis, ether_hydrolysis and nitrile_hydrolysis, +# so it runs by default for all of those families. XY_Addition_MultipleBond is newly +# registered for 'heuristics' and 'crest'; of the two, only 'heuristics' is on by default. ts_adapters = ['heuristics', 'linear', 'AutoTST', 'GCN', 'xtb_gsm', 'orca_neb'] # List here job types to execute by default @@ -587,3 +593,9 @@ def add_rmg_db_candidates(prefix: str) -> None: if path and os.path.isdir(path): RMG_DB_PATH = path break + +# CREST is located once and cached. Standalone builds are looked for under /Local/ce_dana; +# set ARC_CREST_STANDALONE_DIR to another directory to move that search, or to an empty +# string to skip it. +CREST_PATH, CREST_ENV_PATH = get_crest_paths() + diff --git a/arc/species/converter.py b/arc/species/converter.py index 5075ae2279..776f8b9f8e 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -5,7 +5,7 @@ import math import numpy as np import os -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Optional from collections.abc import Iterable from ase import Atoms @@ -50,6 +50,97 @@ ANGL_PRECISION = 0.1 # rad (for both bond angle and dihedral) +def reorder_xyz_string(xyz_str: str, + reverse_atoms: bool = False, + units: str = 'angstrom', + convert_to: str = 'angstrom', + project_directory: Optional[str] = None + ) -> str: + """ + Reorder an XYZ string between ``ATOM X Y Z`` and ``X Y Z ATOM`` with optional unit conversion. + + The order of the input is detected from its first non-empty line, and is preserved when + ``reverse_atoms`` is ``False`` and flipped when it is ``True``. Empty lines are ignored. + A tuple or a list input is joined into a single string with newlines. + + Args: + xyz_str (str): The string xyz format to be converted. + reverse_atoms (bool, optional): Whether to reverse the atoms and coordinates. + units (str, optional): Units of the input coordinates ('angstrom' or 'bohr'). + convert_to (str, optional): The units to convert to (either 'angstrom' or 'bohr'). + project_directory (str, optional): The path to the project directory. + + Raises: + ConverterError: If ``xyz_str`` is not a string, is empty, does not have four space-separated + entries in each of its non-empty lines, has non-float coordinates, + or if ``units`` or ``convert_to`` is neither 'angstrom' nor 'bohr'. + + Returns: str + The converted string xyz format. + """ + if isinstance(xyz_str, tuple): + xyz_str = '\n'.join(xyz_str) + if isinstance(xyz_str, list): + xyz_str = '\n'.join(xyz_str) + if not isinstance(xyz_str, str): + raise ConverterError(f'Expected a string input, got {type(xyz_str)}') + if project_directory is not None: + file_path = os.path.join(project_directory, xyz_str) + if os.path.isfile(file_path): + with open(file_path, 'r') as f: + xyz_str = f.read() + + if units.lower() == 'angstrom' and convert_to.lower() == 'angstrom': + conversion_factor = 1 + elif units.lower() == 'bohr' and convert_to.lower() == 'bohr': + conversion_factor = 1 + elif units.lower() == 'angstrom' and convert_to.lower() == 'bohr': + conversion_factor = constants.angstrom_to_bohr + elif units.lower() == 'bohr' and convert_to.lower() == 'angstrom': + conversion_factor = constants.bohr_to_angstrom + else: + raise ConverterError("Invalid target unit. Choose 'angstrom' or 'bohr'.") + + processed_lines = list() + lxyz = [line for line in xyz_str.strip().splitlines() if line.strip()] + if not lxyz: + raise ConverterError(f'xyz_str has an incorrect format, expected at least one non-empty line, ' + f'got:\n{xyz_str}') + atom_first = not is_str_float(lxyz[0].strip().split()[0]) + + for item in lxyz: + parts = item.strip().split() + + if len(parts) != 4: + raise ConverterError(f'xyz_str has an incorrect format, expected 4 elements in each line, ' + f'got "{item}" in:\n{xyz_str}') + if atom_first: + atom, x_str, y_str, z_str = parts + else: + x_str, y_str, z_str, atom = parts + + try: + x = float(x_str) * conversion_factor + y = float(y_str) * conversion_factor + z = float(z_str) * conversion_factor + + except ValueError as e: + raise ConverterError(f'Could not convert {x_str}, {y_str}, or {z_str} to floats.') from e + + if reverse_atoms and atom_first: + formatted_line = f'{x} {y} {z} {atom}' + elif reverse_atoms and not atom_first: + formatted_line = f'{atom} {x} {y} {z}' + elif not reverse_atoms and atom_first: + formatted_line = f'{atom} {x} {y} {z}' + elif not reverse_atoms and not atom_first: + formatted_line = f'{x} {y} {z} {atom}' + + processed_lines.append(formatted_line) + + return '\n'.join(processed_lines) + + def str_to_xyz(xyz_str: str, project_directory: str | None = None, ) -> dict: diff --git a/arc/species/converter_test.py b/arc/species/converter_test.py index 033bdf78ac..4e9f61ae83 100644 --- a/arc/species/converter_test.py +++ b/arc/species/converter_test.py @@ -7,6 +7,7 @@ import math import os +import shutil import numpy as np from scipy.spatial.transform import Rotation @@ -19,6 +20,7 @@ import arc.species.converter as converter from arc.common import (ARC_PATH, ARC_TESTING_PATH, almost_equal_coords, almost_equal_coords_lists, almost_equal_lists, distance_matrix) +import arc.constants as constants from arc.exceptions import ConverterError from arc.molecule.molecule import Molecule from arc.species.perceive import perceive_molecule_from_xyz @@ -701,6 +703,137 @@ def test_str_to_xyz(self): xyz = converter.str_to_xyz(xyz_format) self.assertEqual(xyz, expected_xyz) + def test_reorder_xyz_string_atom_first(self): + """Test reordering atom-first XYZ strings with unit conversion""" + xyz_format = "C 0.0 1.0 2.0\nH -1.0 0.5 0.0" + converted = converter.reorder_xyz_string(xyz_str=xyz_format, reverse_atoms=True, convert_to="bohr") + converted_lines = converted.splitlines() + self.assertEqual(len(converted_lines), 2) + + x1, y1, z1, s1 = converted_lines[0].split() + self.assertEqual(s1, "C") + self.assertAlmostEqual(float(x1), 0.0) + self.assertAlmostEqual(float(y1), 1.0 * constants.angstrom_to_bohr) + self.assertAlmostEqual(float(z1), 2.0 * constants.angstrom_to_bohr) + + x2, y2, z2, s2 = converted_lines[1].split() + self.assertEqual(s2, "H") + self.assertAlmostEqual(float(x2), -1.0 * constants.angstrom_to_bohr) + self.assertAlmostEqual(float(y2), 0.5 * constants.angstrom_to_bohr) + self.assertAlmostEqual(float(z2), 0.0) + + def test_reorder_xyz_string_coordinate_first(self): + """Test reordering coordinate-first XYZ strings back to atom-last order with conversion""" + xyz_format = "0.0 0.0 0.0 N\n1.0 0.0 0.0 H" + converted = converter.reorder_xyz_string( + xyz_str=xyz_format, + reverse_atoms=False, + units="bohr", + convert_to="angstrom", + ) + # Compared numerically rather than as a byte-exact string: the conversion carries the full + # precision of ``bohr_to_angstrom``, and pinning its digits here would tie a shared physical + # constant to a string literal. + lines = converted.splitlines() + self.assertEqual(len(lines), 2) + self.assertEqual(lines[0], "0.0 0.0 0.0 N") + x, y, z, symbol = lines[1].split() + self.assertEqual(symbol, 'H') + self.assertAlmostEqual(float(x), constants.bohr_to_angstrom, places=8) + self.assertAlmostEqual(float(y), 0.0, places=8) + self.assertAlmostEqual(float(z), 0.0, places=8) + + def test_reorder_xyz_string_all_ordering_branches(self): + """Test all four combinations of input order and reverse_atoms""" + atom_first_str = "C 1.0 2.0 3.0\nH 4.0 5.0 6.0" + coord_first_str = "1.0 2.0 3.0 C\n4.0 5.0 6.0 H" + expected_symbols = ['C', 'H'] + expected_coords = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]] + + atom_first_kept = converter.reorder_xyz_string(xyz_str=atom_first_str, reverse_atoms=False) + coord_first_flipped = converter.reorder_xyz_string(xyz_str=coord_first_str, reverse_atoms=True) + for converted in [atom_first_kept, coord_first_flipped]: + lines = converted.splitlines() + self.assertEqual(len(lines), 2) + for line, symbol, coords in zip(lines, expected_symbols, expected_coords): + tokens = line.split() + self.assertEqual(len(tokens), 4) + self.assertEqual(tokens[0], symbol) + for token, coord in zip(tokens[1:], coords): + self.assertAlmostEqual(float(token), coord, places=8) + + atom_first_flipped = converter.reorder_xyz_string(xyz_str=atom_first_str, reverse_atoms=True) + coord_first_kept = converter.reorder_xyz_string(xyz_str=coord_first_str, reverse_atoms=False) + for converted in [atom_first_flipped, coord_first_kept]: + lines = converted.splitlines() + self.assertEqual(len(lines), 2) + for line, symbol, coords in zip(lines, expected_symbols, expected_coords): + tokens = line.split() + self.assertEqual(len(tokens), 4) + self.assertEqual(tokens[3], symbol) + for token, coord in zip(tokens[:3], coords): + self.assertAlmostEqual(float(token), coord, places=8) + + def test_reorder_xyz_string_identity_unit_conversions(self): + """Test that angstrom to angstrom and bohr to bohr leave the coordinates unchanged""" + xyz_format = "C 1.5 -2.25 3.125\nH -4.0 5.5 0.75" + expected_coords = [[1.5, -2.25, 3.125], [-4.0, 5.5, 0.75]] + for units in ['angstrom', 'bohr']: + converted = converter.reorder_xyz_string(xyz_str=xyz_format, units=units, convert_to=units) + lines = converted.splitlines() + self.assertEqual(len(lines), 2) + for line, coords in zip(lines, expected_coords): + tokens = line.split() + for token, coord in zip(tokens[1:], coords): + self.assertAlmostEqual(float(token), coord, places=8) + + def test_reorder_xyz_string_errors_and_contracts(self): + """Test the documented error cases, the blank line contract, and tuple/list inputs""" + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str="C 0.0 0.0 0.0", convert_to="furlong") + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str="C 0.0 0.0 0.0", units="furlong") + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str="C 0.0 0.0 0.0\nH 0.0 0.0") + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str="C 0.0 0.0 spam") + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str=5) + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str={'symbols': ('C',)}) + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str='') + with self.assertRaises(ConverterError): + converter.reorder_xyz_string(xyz_str='\n \n') + + blank_line_input = 'C 0.0 0.0 0.0\n\nH 0.0 0.0 1.0' + self.assertEqual(converter.reorder_xyz_string(xyz_str=blank_line_input).splitlines(), + ['C 0.0 0.0 0.0', 'H 0.0 0.0 1.0']) + + joined_from_tuple = converter.reorder_xyz_string(xyz_str=('C 0.0 0.0 0.0', 'H 0.0 0.0 1.0')) + joined_from_list = converter.reorder_xyz_string(xyz_str=['C 0.0 0.0 0.0', 'H 0.0 0.0 1.0']) + self.assertEqual(joined_from_tuple, 'C 0.0 0.0 0.0\nH 0.0 0.0 1.0') + self.assertEqual(joined_from_list, joined_from_tuple) + + def test_reorder_xyz_string_reads_from_project_directory(self): + """Test that a file name relative to the project directory is read, and a plain xyz string is not""" + project_directory = os.path.join(ARC_PATH, 'Projects', 'arc_project_for_testing_delete_after_usage_reorder') + self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True) + os.makedirs(project_directory, exist_ok=True) + file_name = 'coords.xyz' + with open(os.path.join(project_directory, file_name), 'w') as f: + f.write('C 1.0 2.0 3.0\nH 4.0 5.0 6.0\n') + self.assertEqual(converter.reorder_xyz_string(xyz_str=file_name, project_directory=project_directory), + 'C 1.0 2.0 3.0\nH 4.0 5.0 6.0') + self.assertEqual(converter.reorder_xyz_string(xyz_str='O 0.0 0.0 0.0', project_directory=project_directory), + 'O 0.0 0.0 0.0') + + def test_angstrom_bohr_conversion_constants(self): + """Test the angstrom/bohr conversion constants against known values, independently of the code that uses them""" + self.assertAlmostEqual(constants.angstrom_to_bohr, 1.8897261246, places=8) + self.assertAlmostEqual(constants.bohr_to_angstrom, 0.5291772109, places=8) + self.assertAlmostEqual(constants.angstrom_to_bohr * constants.bohr_to_angstrom, 1.0, places=12) + def test_xyz_to_str(self): """Test converting an ARC xyz format to a string xyz format""" xyz_str1 = converter.xyz_to_str(xyz_dict=self.xyz1['dict']) diff --git a/arc/species/species_test.py b/arc/species/species_test.py index 1abac6ca85..1fa07fedc3 100644 --- a/arc/species/species_test.py +++ b/arc/species/species_test.py @@ -3614,6 +3614,7 @@ def test_from_dict(self): ts_dict = self.tsg1.as_dict() tsg = TSGuess(ts_dict=ts_dict) self.assertEqual(tsg.method, 'autotst') + self.assertEqual(tsg.method_sources, ['autotst']) ts_dict_for_report = self.tsg1.as_dict(for_report=True) self.assertEqual(list(ts_dict_for_report.keys()), ['method', 'method_sources', 'method_index', 'success', 'index', 'conformer_index', 'initial_xyz', 'opt_xyz']) diff --git a/devtools/crest_environment.yml b/devtools/crest_environment.yml new file mode 100644 index 0000000000..2291e72d37 --- /dev/null +++ b/devtools/crest_environment.yml @@ -0,0 +1,6 @@ +name: crest_env +channels: + - conda-forge +dependencies: + - python>=3.7 + - crest=2.12 diff --git a/devtools/install_crest.sh b/devtools/install_crest.sh new file mode 100644 index 0000000000..1086ec9db2 --- /dev/null +++ b/devtools/install_crest.sh @@ -0,0 +1,64 @@ +#!/bin/bash -l +set -eo pipefail + +if command -v micromamba &> /dev/null; then + echo "✔️ Micromamba is installed." + COMMAND_PKG=micromamba +elif command -v mamba &> /dev/null; then + echo "✔️ Mamba is installed." + COMMAND_PKG=mamba +elif command -v conda &> /dev/null; then + echo "✔️ Conda is installed." + COMMAND_PKG=conda +else + echo "❌ Micromamba, Mamba, or Conda is required. Please install one." + exit 1 +fi + +if [ "$COMMAND_PKG" = "micromamba" ]; then + eval "$(micromamba shell hook --shell=bash)" +else + BASE=$(conda info --base) + . "$BASE/etc/profile.d/conda.sh" +fi + +ENV_FILE="devtools/crest_environment.yml" + +if [ ! -f "$ENV_FILE" ]; then + echo "❌ File not found: $ENV_FILE" + exit 1 +fi + +if $COMMAND_PKG env list | grep -q '^crest_env\s'; then + echo ">>> Updating existing crest_env..." + $COMMAND_PKG env update -n crest_env -f "$ENV_FILE" --prune +else + echo ">>> Creating new crest_env..." + $COMMAND_PKG env create -n crest_env -f "$ENV_FILE" -y +fi + +echo ">>> Checking CREST installation..." + +if [ "$COMMAND_PKG" = "micromamba" ]; then + CREST_RUNNER="micromamba run -n crest_env" + CREST_LISTER="micromamba list -n crest_env" +else + CREST_RUNNER="conda run -n crest_env" + CREST_LISTER="conda list -n crest_env" +fi + +if $CREST_RUNNER crest --version &> /dev/null; then + version_output=$($CREST_RUNNER crest --version 2>&1) + echo "$version_output" + installed_version=$(printf '%s' "$version_output" | tr '\n' ' ' | sed -n 's/.*Version[[:space:]]\+\([0-9.][0-9.]*\).*/\1/p') + if [ "$installed_version" != "2.12" ]; then + echo "❌ CREST version mismatch (expected 2.12)." + exit 1 + fi + echo "✔️ CREST 2.12 is successfully installed." +else + echo "❌ CREST is not found in PATH. Please check the environment." + exit 1 +fi + +echo "✅ Done installing CREST (crest_env)." diff --git a/docs/source/TS_search.rst b/docs/source/TS_search.rst index a2995372fa..8b0028cf51 100644 --- a/docs/source/TS_search.rst +++ b/docs/source/TS_search.rst @@ -493,4 +493,125 @@ DOI `10.26434/chemrxiv.15001681/v1 `_. +CREST +^^^^^ + +CREST is an external conformational sampling tool used by ARC as a TS-search wrapper stage. +In ARC's current flow, CREST is applied to TS seeds generated by base TS search methods and uses +family-specific constraints from ARC. + +CREST is **opt-in only** — it needs its own environment and an external binary, so it is +intentionally absent from the default ``ts_adapters`` list. Install it once with +``make install-crest`` and request it per run: + +.. code-block:: yaml + + ts_adapters: + - heuristics + - crest + +Current ARC family support for CREST: + +- ``H_Abstraction`` (RMG family reference: + `H_Abstraction `_). +- ``XY_Addition_MultipleBond``, seeded by the four-centre XY-addition builder. This family belongs + to RMG's ``halogens`` family set rather than the ``default`` set, so it is only classified — and + the builder only reached — when that set is enabled for the run, i.e. with + ``rmg_family_set: 'halogens'`` or ``rmg_family_set: 'all'``. A list of family-set labels such as + ``['default', 'halogens']`` does not work: ``get_all_families()`` returns such a list unexpanded, + so the labels are taken for family names and every RMG family is dropped. +- ``carbonyl_based_hydrolysis``, ``ether_hydrolysis`` and ``nitrile_hydrolysis``, constrained on all + six pairwise distances among the four reactive atoms (electrophilic centre, leaving group, water + oxygen and the transferring water hydrogen). Four atoms have six internal degrees of freedom, so + a smaller set would leave the reactive core free to relax away from the saddle. These are + ARC-native families and are always available. + +Each TS seed becomes its own CREST metadynamics job, and a heuristics sweep can produce dozens of +seeds for a single reaction. At most ``MAX_CREST_SEEDS`` of the most geometrically distinct seeds +are submitted; the number dropped is logged. + +Locating the binary: standalone builds are searched for under ``/Local/ce_dana``. Set +``ARC_CREST_STANDALONE_DIR`` to search a different directory, or to an empty string to skip that +search and rely on a conda/mamba ``crest_env`` or ``PATH``. + +Upgrading an existing checkout: ``arc/constants.pxd`` gained ``angstrom_to_bohr``, so a stale +compiled ``arc.constants`` shadows the source and lacks the symbol. Re-run ``make compile`` after +pulling, otherwise the first CREST job fails with an ``AttributeError`` rather than at import. + +External references: + +- `CREST documentation `_ +- `CREST constrained sampling example `_ + +Wrapper Extension Guide +""""""""""""""""""""""" + +Use this guide when extending CREST-based TS workflows in ARC (for example, adding hydrolysis support to CREST, +or allowing CREST to wrap a new TS seed source adapter). + +ARC uses a neutral wrapper hub API for TS seed generation and wrapper-specific constraints: + +- ``arc.job.adapters.ts.seed_hub.get_ts_seeds(...)`` +- ``arc.job.adapters.ts.seed_hub.get_wrapper_constraints(...)`` + +Current status +"""""""""""""" + +- ``CrestAdapter`` requests seeds using ``base_adapter='heuristics'``. +- ``CrestAdapter`` requests constraints using ``wrapper='crest'``. +- CREST constraints are implemented for ``H_Abstraction``, ``XY_Addition_MultipleBond`` and the + three hydrolysis families. +- Constraints pin interatomic distances only. An angle restraint is redundant once the distances + spanning the reactive core are pinned, and is ill-conditioned near linearity. + +Seed schema contract +"""""""""""""""""""" + +``get_ts_seeds(...)`` returns a list of seed dictionaries with the following fields: + +- ``xyz``: Cartesian coordinates dictionary. +- ``family``: Reaction family associated with the seed. +- ``method``: Method label for provenance. +- ``source_adapter``: TS-search adapter id that generated the seed. +- ``metadata``: Optional adapter-specific metadata dictionary. + +Extension instructions: Add a new family to CREST +""""""""""""""""""""""""""""""""""""""""""""""""" + +1. Update ``get_ts_seeds(...)`` logic in ``arc/job/adapters/ts/seed_hub.py`` only if the seed generation path changes. +2. Add family-specific CREST constraints in ``_get_crest_constraints(...)`` (or family helper it calls) in + ``arc/job/adapters/ts/seed_hub.py``. +3. Add/update tests in ``arc/job/adapters/ts/heuristics_test.py`` (``TestHeuristicsHub``). +4. Update ``ts_adapters_by_rmg_family`` mapping if CREST should be enabled for that family. + +Extension instructions: Let CREST wrap a new TS seed adapter +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +1. Add a ``base_adapter`` branch in ``get_ts_seeds(...)``. +2. Ensure the returned seed objects satisfy the seed schema contract. +3. Reuse ``get_wrapper_constraints(wrapper='crest', ...)`` with those seeds. +4. Add tests for the new adapter branch and constraints compatibility. + +Minimal usage pattern +""""""""""""""""""""" + +.. code-block:: python + + from arc.job.adapters.ts.seed_hub import get_ts_seeds, get_wrapper_constraints + + seeds = get_ts_seeds( + reaction=rxn, + base_adapter='heuristics', + dihedral_increment=30, + ) + for seed in seeds: + crest_constraints = get_wrapper_constraints( + wrapper='crest', + reaction=rxn, + seed=seed, + ) + if crest_constraints is None: + continue + # run CREST with crest_constraints["A"], crest_constraints["H"], crest_constraints["B"] + .. include:: links.txt