diff --git a/arc/job/local.py b/arc/job/local.py index c4a2928215..7eb50bf388 100644 --- a/arc/job/local.py +++ b/arc/job/local.py @@ -8,6 +8,7 @@ import os import re import shutil +import signal import subprocess import time @@ -28,6 +29,7 @@ def execute_command(command: str | list[str], shell: bool = True, no_fail: bool = False, executable: str | None = None, + timeout: float | None = None, ) -> tuple[list | None, list | None]: """ Execute a command. @@ -35,6 +37,8 @@ def execute_command(command: str | list[str], Notes: If ``no_fail`` is ``True``, then a warning is logged and ``False`` is returned so that the calling function can debug the situation. + A command that exceeds ``timeout`` is not retried: the point of a deadline is to bound + the call, so the timeout is reported through the same channels as a terminal failure. Args: command (str | list[str]): An array of string commands to send. @@ -42,6 +46,20 @@ def execute_command(command: str | list[str], no_fail (bool, optional): If ``True`` then ARC will not crash if an error is encountered. executable (str, optional): Select a specific shell to run with, e.g., '/bin/bash'. Default shell of the subprocess command is '/bin/sh'. + timeout (float, optional): The number of seconds to wait for the command to complete. + If the command does not complete in time, it is killed along with + every process it spawned that stayed in its process group. A + descendant that puts itself in a group of its own, e.g. via + ``setsid``, is outside what signalling a group can reach and + survives. ``None``, the default, waits forever. + Note that the deadline is not exact: a timing out call returns + after ``timeout`` plus up to twice the grace period of + ``_kill_process_group()``, i.e. up to 10 seconds late by default. + + Raises: + SettingsError: If the command timed out and ``no_fail`` is ``False``. A non-zero exit status + is not an error here: the command is run without ``check=True``, so a failing + command's stderr is returned as output rather than raised. Returns: tuple[list, list]: - A list of lines of standard output stream. @@ -55,12 +73,28 @@ def execute_command(command: str | list[str], sleep_time = 60 # Seconds while i < max_times_to_try: try: + if timeout is not None: + stdout, stderr = _run_with_timeout(command=command, shell=shell, + executable=executable, timeout=timeout) + return _format_stdout(stdout), _format_stdout(stderr) if executable is None: completed_process = subprocess.run(command, shell=shell, capture_output=True) else: completed_process = subprocess.run(command, shell=shell, capture_output=True, executable=executable) return _format_stdout(completed_process.stdout), _format_stdout(completed_process.stderr) + except subprocess.TimeoutExpired: + message = f'The command "{command}" did not complete within {timeout} seconds ' \ + f'and was terminated along with all of the processes it spawned.' + if no_fail: + logger.warning(message) + return None, None + logger.error(message) + raise SettingsError(f'{message}\nConsider increasing the timeout, or check whether this is a ' + f'server issue by executing the command manually on the server.') from None except subprocess.CalledProcessError as e: + # Note: ``subprocess.run()`` is called without ``check=True`` and therefore never + # raises this, so this handler and the retry loop it drives are dead code. They are + # left as they are; the ``TimeoutExpired`` handler above is the only live one. error = e # Store the error so we can raise a SettingsError if needed. if no_fail: _output_command_error_message(command, e, logger.warning) @@ -82,6 +116,110 @@ def execute_command(command: str | list[str], f'sbatch path required in the submit_command dictionary.') +def _run_with_timeout(command: list[str], + shell: bool, + executable: str | None, + timeout: float, + ) -> tuple[bytes, bytes]: + """ + Run a command under a deadline, killing its process group if the deadline expires. + + The command is started in its own session (hence its own process group) so that the whole tree + can be signalled at once. Killing only the direct child is not enough: with ``shell=True`` the + direct child is a shell that commonly ``exec``s its last command, and anything it backgrounded + would survive as an orphan. + + The containment this buys reaches exactly as far as the process group does. A descendant that + leaves the group, e.g. by calling ``setsid`` itself, receives neither signal and keeps running; + holding it too would take a cgroup, which is beyond what this function attempts. + + Note that a timing out call returns only after the grace periods of ``_kill_process_group()``, + which sleeps for one before escalating to SIGKILL and then waits up to another for the child. + The call is therefore bounded by ``timeout + 2 * grace_period``, i.e. up to 10 seconds beyond + the deadline with the default grace period, rather than by ``timeout`` alone. + + Args: + command (list[str]): The command to run, already normalized by ``execute_command()`` into + the single-element list that ``subprocess`` expects. + shell (bool): Specifies whether the command should be executed using bash instead of Python. + executable (str | None): Select a specific shell to run with, e.g., '/bin/bash'. + timeout (float): The number of seconds to wait for the command to complete. + + Raises: + subprocess.TimeoutExpired: If the command did not complete within ``timeout`` seconds. + + Returns: tuple[bytes, bytes]: + - The standard output stream. + - The standard error stream. + """ + process = subprocess.Popen(command, + shell=shell, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + executable=executable, + start_new_session=True, + ) + try: + return process.communicate(timeout=timeout) + except subprocess.TimeoutExpired: + _kill_process_group(process=process) + raise + finally: + # ``Popen`` is deliberately not used as a context manager here: its ``__exit__()`` waits for + # the child without a timeout, which would defeat the very deadline this function enforces + # whenever a child cannot be reaped promptly. Closing the pipes is therefore done here, and + # reaping is left to ``communicate()`` on the normal path and to the bounded wait of + # ``_kill_process_group()`` on the timeout path. + for stream in (process.stdout, process.stderr): + if stream is not None: + stream.close() + + +def _kill_process_group(process: subprocess.Popen, + grace_period: float = 5, + ) -> None: + """ + Kill a process and every process it spawned that is still in its process group. + + SIGTERM is sent first to let the tree shut down cleanly, then SIGKILL after ``grace_period``. + Falls back to killing just the given process if it does not have a process group of its own, + which also guards against ever signalling the group ARC itself runs in. + + Args: + process (subprocess.Popen): The process to kill, started with ``start_new_session=True``. + grace_period (float, optional): The number of seconds to allow the tree to exit on SIGTERM. + """ + try: + pgid = os.getpgid(process.pid) + except OSError: + pgid = None + if pgid is None or pgid == os.getpgid(0): + process.kill() + else: + try: + os.killpg(pgid, signal.SIGTERM) + except OSError: + pass # The group is already gone, which is the outcome this call is after anyway. + # The child is deliberately not reaped here. It is the group leader, so while it remains + # unreaped its pid cannot be recycled and ``pgid`` stays pinned to this group, which is + # what makes the SIGKILL below safe to send. Reaping first would open a window in which + # the pid is free and the SIGKILL could land on an unrelated process group. + # The escalation also cannot be conditioned on the child having survived SIGTERM: with + # ``shell=True`` the direct child is a shell that dies on SIGTERM, while a grandchild + # that ignores SIGTERM is precisely the process that gets orphaned. + time.sleep(grace_period) + try: + os.killpg(pgid, signal.SIGKILL) + except OSError: + pass # The group exited on SIGTERM, which is the preferred outcome. + try: + process.wait(timeout=grace_period) + except subprocess.TimeoutExpired: + logger.warning(f'Could not terminate process {process.pid} after it timed out. ' + f'It is left unreaped rather than waited for indefinitely, so that a child ' + f'which cannot be killed does not also hang ARC.') + + def _output_command_error_message(command: list[str], error: subprocess.CalledProcessError, logging_func, diff --git a/arc/job/local_test.py b/arc/job/local_test.py index 68d4bd496f..9f89e4e14a 100644 --- a/arc/job/local_test.py +++ b/arc/job/local_test.py @@ -8,11 +8,71 @@ import datetime import os import shutil +import signal +import sys +import tempfile +import time import unittest from unittest.mock import patch import arc.job.local as local from arc.common import ARC_PATH +from arc.exceptions import SettingsError + + +def process_is_alive(pid: int) -> bool: + """ + Whether a process with the given pid exists and is not a zombie. + + Args: + pid (int): The process ID to check. + + Returns: + bool: Whether the process is alive. + """ + try: + os.kill(pid, 0) + except OSError: + return False + try: + with open(f'/proc/{pid}/stat', 'r') as f: + return f.read().rsplit(')', 1)[-1].split()[0] != 'Z' + except (OSError, IndexError): + return True + + +def get_parent_pid(pid: int) -> int | None: + """ + Get the parent pid of a process, used to report whether a process was reparented to init. + + Args: + pid (int): The process ID to check. + + Returns: + int | None: The parent process ID, or ``None`` if it could not be determined. + """ + try: + with open(f'/proc/{pid}/stat', 'r') as f: + return int(f.read().rsplit(')', 1)[-1].split()[1]) + except (OSError, IndexError, ValueError): + return None + + +def _read_pid(path: str) -> int | None: + """ + Read a pid from a file, returning ``None`` if it is absent or unreadable. + + Args: + path (str): The path of the file holding the pid. + + Returns: + int | None: The pid, or ``None``. + """ + try: + with open(path, 'r') as f: + return int(f.read().strip()) + except (OSError, ValueError): + return None class TestLocal(unittest.TestCase): @@ -37,6 +97,72 @@ def test_execute_command(self): self.assertIn('adapter.py', out1[0]) self.assertIn('ssh.py', out1[0]) + def test_execute_command_without_a_timeout_is_unchanged(self): + """Test that not passing a timeout leaves the subprocess call exactly as it was""" + with patch('arc.job.local.subprocess.run') as mock_run: + mock_run.return_value.stdout, mock_run.return_value.stderr = b'ok\n', b'' + local.execute_command('ls') + local.execute_command('ls', executable='/bin/bash') + self.assertEqual(len(mock_run.call_args_list), 2) + for call in mock_run.call_args_list: + self.assertNotIn('timeout', call.kwargs) + self.assertNotIn('start_new_session', call.kwargs) + self.assertEqual(mock_run.call_args_list[0].kwargs, {'shell': True, 'capture_output': True}) + self.assertEqual(mock_run.call_args_list[1].kwargs, + {'shell': True, 'capture_output': True, 'executable': '/bin/bash'}) + + def test_execute_command_with_a_generous_timeout(self): + """Test that a command which completes in time is unaffected by a timeout""" + stdout, stderr = local.execute_command('echo hello', timeout=60) + self.assertEqual(stdout, ['hello']) + self.assertEqual(stderr, []) + + def test_execute_command_timeout_kills_spawned_children(self): + """Test that a timing out command is killed along with the processes it spawned""" + # The grandchild ignores SIGTERM, so this pins both halves of the kill: + # signalling the whole process group rather than only the direct child, and escalating + # to SIGKILL. A grandchild that dies on SIGTERM cannot tell the escalation apart. + # The direct child here is a plain ``sleep`` that does die on SIGTERM, which is exactly + # why the escalation must not be conditioned on the direct child having survived. + temp_dir = tempfile.mkdtemp() + pid_path = os.path.join(temp_dir, 'grandchild.pid') + script_path = os.path.join(temp_dir, 'stubborn_grandchild.py') + with open(script_path, 'w') as f: + f.write('import os, signal, sys, time\n' + 'signal.signal(signal.SIGTERM, signal.SIG_IGN)\n' + "with open(sys.argv[1], 'w') as f:\n" + " f.write(str(os.getpid()))\n" + 'time.sleep(300)\n') + # The shell waits for the grandchild to record its pid before starting the long sleep, so a + # loaded worker cannot have the timeout fire while the grandchild is still starting up and + # fail the assertion below even though process group cleanup worked correctly. + command = f'{sys.executable} {script_path} {pid_path} & ' \ + f'while [ ! -s {pid_path} ]; do sleep 0.05; done; sleep 300' + try: + with self.assertRaises(SettingsError): + local.execute_command(command, timeout=3) + self.assertTrue(os.path.isfile(pid_path), 'The grandchild never recorded its pid.') + with open(pid_path, 'r') as f: + grandchild_pid = int(f.read().strip()) + for _ in range(200): + if not process_is_alive(grandchild_pid): + break + time.sleep(0.1) + self.assertFalse(process_is_alive(grandchild_pid), + f'The spawned process {grandchild_pid} was orphaned instead of killed ' + f'(parent pid {get_parent_pid(grandchild_pid)}).') + finally: + for pid in [_read_pid(pid_path)]: + if pid is not None and process_is_alive(pid): + os.kill(pid, signal.SIGKILL) + shutil.rmtree(temp_dir, ignore_errors=True) + + def test_execute_command_timeout_with_no_fail(self): + """Test that a timing out command returns None, None when no_fail is True""" + stdout, stderr = local.execute_command('sleep 300', no_fail=True, timeout=2) + self.assertIsNone(stdout) + self.assertIsNone(stderr) + def test_determine_job_id(self): """Test determining a job ID from the stdout of a job submission command.""" # Slurm diff --git a/arc/reaction/reaction.py b/arc/reaction/reaction.py index f6a411a5dc..21a518eb93 100644 --- a/arc/reaction/reaction.py +++ b/arc/reaction/reaction.py @@ -154,9 +154,18 @@ def __init__(self, @property def atom_map(self): - """The reactants to products atom map""" + """ + The reactants to products atom map. + + Warning: + Mapping requires 3D coordinates, so reading this property generates a cheap conformer + (an RDKit embedding followed by a force field optimization) for every reactant and + product that has none, storing it on the species. Until check_atom_balance() stopped + generating a conformer for every species this was guaranteed to be a no-op here, and + the cost merely moved rather than appeared. + """ if self._atom_map is None \ - and all(species.get_xyz(generate=False) is not None for species in self.r_species + self.p_species): + and all(species.get_xyz(generate=True) is not None for species in self.r_species + self.p_species): _atom_map = map_reaction(rxn=self, backend='ARC') if _atom_map is not None: self._atom_map = _atom_map @@ -770,18 +779,18 @@ def check_atom_balance(self, for reactant in self.r_species: count = self.get_species_count(species=reactant, well=0) - xyz = reactant.get_xyz(generate=True) - if xyz is not None and xyz: - r_well += (xyz_to_str(xyz) + '\n') * count + entry = _get_atom_balance_entry(species=reactant) + if entry: + r_well += (entry + '\n') * count else: r_well = '' break for product in self.p_species: count = self.get_species_count(species=product, well=1) - xyz = product.get_xyz(generate=True) - if xyz is not None and xyz: - p_well += (xyz_to_str(xyz) + '\n') * count + entry = _get_atom_balance_entry(species=product) + if entry: + p_well += (entry + '\n') * count else: p_well = '' break @@ -1258,6 +1267,47 @@ def get_rxn_smiles(self) -> str | None: return ".".join(smiles_r)+">>"+".".join(smiles_p) +def _get_atom_balance_entry(species: ARCSpecies) -> str: + """ + Get an xyz string representation of a species to be used in an atom balance check. + + An atom balance check only counts element symbols, it never uses the coordinates. + Therefore, an already available geometry is used if the species has one, and otherwise the + element symbols are read off the species' 2D graph and given placeholder coordinates. + This avoids cheaply generating a 3D conformer (an RDKit embedding followed by a force field + optimization, by far the most expensive step of an atom balance check) merely to count atoms. + + Warning: + In the 2D graph case every atom is placed at the origin, so the returned string is a + fabricated geometry. It must only ever be used for counting elements, which is why this + function is private to this module. + + Note: + The conditions under which coordinates and a 2D graph are considered mirror + ``ARCSpecies.get_xyz(generate=True)``: a TS species only ever reports the coordinates of + its TS guesses, the ``mol_list`` fallback matches that method's + ``self.mol is not None or self.mol_list is not None`` guard, and the resulting element + counts are identical to those of a generated conformer, which is embedded from the very + same 2D graph. + + Args: + species (ARCSpecies): The species to represent. + + Returns: + str: An xyz string representation of the species, + or an empty string if the species has neither coordinates nor a 2D graph. + """ + xyz = species.get_xyz(generate=False) + if xyz is not None and xyz: + return xyz_to_str(xyz) + if species.is_ts: + return '' + mol = species.mol if species.mol is not None else (species.mol_list[0] if species.mol_list else None) + if mol is None or not len(mol.atoms): + return '' + return '\n'.join(f'{atom.element.symbol} 0.0 0.0 0.0' for atom in mol.atoms) + + def remove_dup_species(species_list: list[ARCSpecies]) -> list[ARCSpecies]: """ Remove duplicate species from a species list. diff --git a/arc/reaction/reaction_test.py b/arc/reaction/reaction_test.py index 9870e6b80e..6d31baa6e8 100644 --- a/arc/reaction/reaction_test.py +++ b/arc/reaction/reaction_test.py @@ -17,9 +17,10 @@ from arc.family.family import get_all_families, get_rmg_recommended_family_sets from arc.imports import settings from arc.main import ARC -from arc.reaction.reaction import ARCReaction, remove_dup_species +from arc.reaction.reaction import ARCReaction, _get_atom_balance_entry, remove_dup_species from arc.scheduler import Scheduler from arc.species import ARCSpecies +from arc.species.converter import xyz_to_str from arc.mapping.engine import check_atom_map @@ -191,10 +192,6 @@ def test_as_dict(self): 'label': 'CH4 + OH <=> CH3 + H2O', 'multiplicity': 2, 'p_species': [{'bond_corrections': {'C-H': 3}, - 'cheap_conformer': 'C 0.00000000 0.00000001 -0.00000000\n' - 'H 1.06690511 -0.17519582 0.05416493\n' - 'H -0.68531716 -0.83753536 -0.02808565\n' - 'H -0.38158795 1.01273118 -0.02607927', 'label': 'CH3', 'long_thermo_description': "Bond corrections: {'C-H': 3}\n", 'mol': {'atom_order': rxn_dict_1['p_species'][0]['mol']['atom_order'], @@ -204,9 +201,6 @@ def test_as_dict(self): 'multiplicity': 2, 'number_of_rotors': 0}, {'bond_corrections': {'H-O': 2}, - 'cheap_conformer': 'O -0.00032832 0.39781490 0.00000000\n' - 'H -0.76330345 -0.19953755 0.00000000\n' - 'H 0.76363177 -0.19827735 0.00000000', 'label': 'H2O', 'long_thermo_description': "Bond corrections: {'H-O': 2}\n", 'mol': {'atom_order': rxn_dict_1['p_species'][1]['mol']['atom_order'], @@ -217,11 +211,6 @@ def test_as_dict(self): 'number_of_rotors': 0}], 'products': ['CH3', 'H2O'], 'r_species': [{'bond_corrections': {'C-H': 4}, - 'cheap_conformer': 'C -0.00000000 -0.00000000 0.00000000\n' - 'H -0.63306457 -0.78034118 -0.42801448\n' - 'H -0.38919244 0.98049560 -0.28294367\n' - 'H 0.00329661 -0.09013273 1.08846898\n' - 'H 1.01896040 -0.11002169 -0.37751083', 'label': 'CH4', 'long_thermo_description': "Bond corrections: {'C-H': 4}\n", 'mol': {'atom_order': rxn_dict_1['r_species'][0]['mol']['atom_order'], @@ -231,8 +220,6 @@ def test_as_dict(self): 'multiplicity': 1, 'number_of_rotors': 0}, {'bond_corrections': {'H-O': 1}, - 'cheap_conformer': 'O 0.00000000 0.00000000 0.61310000\n' - 'H 0.00000000 0.00000000 -0.61310000', 'label': 'OH', 'long_thermo_description': "Bond corrections: {'H-O': 1}\n", 'mol': {'atom_order': rxn_dict_1['r_species'][1]['mol']['atom_order'], @@ -251,10 +238,6 @@ def test_as_dict(self): 'family': 'Disproportionation', 'multiplicity': 1, 'p_species': [{'bond_corrections': {'H-N': 3}, - 'cheap_conformer': 'N 0.00064924 -0.00099698 0.29559292\n' - 'H -0.41786606 0.84210396 -0.09477452\n' - 'H -0.52039228 -0.78225292 -0.10002797\n' - 'H 0.93760911 -0.05885406 -0.10079043', 'label': 'NH3', 'long_thermo_description': "Bond corrections: {'H-N': 3}\n", 'mol': {'atom_order': rxn_dict_6['p_species'][0]['mol']['atom_order'], @@ -269,10 +252,6 @@ def test_as_dict(self): 3 H u0 p0 c0 {1,S} 4 H u0 p0 c0 {1,S}""", 'bond_corrections': {'H-N': 2, 'N=N': 1}, - 'cheap_conformer': 'N -0.09608641 0.00717098 -0.00429305\n' - 'N 1.31984473 -0.09850040 -0.31487335\n' - 'H -0.59122841 -0.74658751 0.47254546\n' - 'H -0.63252990 0.83791693 -0.25485633', 'label': 'H2NN[S]', 'long_thermo_description': rxn_dict_6['p_species'][1]['long_thermo_description'], 'mol': {'atom_order': rxn_dict_6['p_species'][1]['mol']['atom_order'], @@ -284,9 +263,6 @@ def test_as_dict(self): 'original_label': 'H2NN(S)'}], 'products': ['H2NN[S]', 'NH3'], 'r_species': [{'bond_corrections': {'H-N': 2}, - 'cheap_conformer': 'N 0.00016375 0.40059499 0.00000000\n' - 'H -0.83170922 -0.19995756 0.00000000\n' - 'H 0.83154548 -0.20063742 0.00000000', 'label': 'NH2', 'long_thermo_description': "Bond corrections: {'H-N': 2}\n", 'mol': {'atom_order': rxn_dict_6['r_species'][0]['mol']['atom_order'], @@ -296,11 +272,6 @@ def test_as_dict(self): 'multiplicity': 2, 'number_of_rotors': 0}, {'bond_corrections': {'H-N': 3, 'N-N': 1}, - 'cheap_conformer': 'N -0.46751749 0.03795671 0.31180026\n' - 'N 0.79325823 -0.46038094 -0.24114357\n' - 'H -1.19307188 -0.63034971 0.05027053\n' - 'H -0.69753009 0.90231202 -0.17907452\n' - 'H 1.56486123 0.15046192 0.05814730', 'label': 'N2H3', 'long_thermo_description': rxn_dict_6['r_species'][1]['long_thermo_description'], 'mol': {'atom_order': rxn_dict_6['r_species'][1]['mol']['atom_order'], @@ -730,6 +701,206 @@ def test_check_atom_balance(self): p_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='H2O', smiles='O')]) + def test_check_atom_balance_does_not_generate_conformers(self): + """Test that the Reaction check_atom_balance method does not generate a 3D conformer""" + def explode(*args, **kwargs): + raise AssertionError('check_atom_balance generated a 3D conformer') + + original_get_cheap_conformer = ARCSpecies.get_cheap_conformer + original_generate_conformers = ARCSpecies.generate_conformers + ARCSpecies.get_cheap_conformer, ARCSpecies.generate_conformers = explode, explode + try: + ch4 = ARCSpecies(label='CH4', smiles='C') + oh = ARCSpecies(label='OH', smiles='[OH]') + ch3 = ARCSpecies(label='CH3', smiles='[CH3]') + h2o = ARCSpecies(label='H2O', smiles='O') + rxn = ARCReaction(r_species=[ch4, oh], p_species=[ch3, h2o]) + self.assertTrue(rxn.check_atom_balance()) + # None of the species were given coordinates, and none were generated either. + for spc in [ch4, oh, ch3, h2o]: + self.assertIsNone(spc.cheap_conformer) + self.assertIsNone(spc.get_xyz(generate=False)) + # An imbalance is still detected without any coordinates. + with self.assertRaises(ReactionError): + ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='H2O', smiles='O')]) + finally: + ARCSpecies.get_cheap_conformer = original_get_cheap_conformer + ARCSpecies.generate_conformers = original_generate_conformers + + def test_check_atom_balance_when_conformer_generation_fails(self): + """Test that an imbalance is caught even for a species no conformer can be generated for. + + A species whose graph is valid but which the force field cannot embed, e.g. a strained cage, + leaves get_xyz(generate=True) returning None. The check used to empty the well and pass such + a reaction silently; it now reads the element symbols off the 2D graph, which is what they + were going to be counted from anyway, and catches the imbalance. + """ + def fail_to_generate(*args, **kwargs): + return None + + original_get_cheap_conformer = ARCSpecies.get_cheap_conformer + original_generate_conformers = ARCSpecies.generate_conformers + ARCSpecies.get_cheap_conformer, ARCSpecies.generate_conformers = fail_to_generate, fail_to_generate + try: + unembeddable = ARCSpecies(label='CH4', smiles='C') + self.assertIsNone(unembeddable.get_xyz(generate=True)) + with self.assertRaises(ReactionError): + ARCReaction(r_species=[unembeddable, ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='H2O', smiles='O')]) + # A balanced reaction built from the same unembeddable species is still accepted. + rxn = 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')]) + self.assertTrue(rxn.check_atom_balance()) + finally: + ARCSpecies.get_cheap_conformer = original_get_cheap_conformer + ARCSpecies.generate_conformers = original_generate_conformers + + def test_check_atom_balance_species_shapes(self): + """Test the Reaction check_atom_balance method for all species shapes that reach it""" + h2o_xyz = {'symbols': ('O', 'H', 'H'), 'isotopes': (16, 1, 1), + 'coords': ((0.0, 0.0, 0.1), (0.0, 0.8, -0.5), (0.0, -0.8, -0.5))} + ch4_adjlist = """1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S}""" + + # 1. Species with a 2D graph and no coordinates (the common case, e.g., from SMILES). + rxn = 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')]) + self.assertTrue(rxn.check_atom_balance()) + + # 2. Species defined by an adjacency list. + rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', adjlist=ch4_adjlist), + ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]'), ARCSpecies(label='H2O', smiles='O')]) + self.assertTrue(rxn.check_atom_balance()) + + # 3. A species defined only by coordinates, with no 2D graph at all. + h2o_no_mol = ARCSpecies(label='H2O', xyz=h2o_xyz) + h2o_no_mol.mol, h2o_no_mol.mol_list = None, None + self.assertIsNone(h2o_no_mol.mol) + rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]'), h2o_no_mol]) + self.assertTrue(rxn.check_atom_balance()) + h2o_no_mol_2 = ARCSpecies(label='H2O', xyz=h2o_xyz) + h2o_no_mol_2.mol, h2o_no_mol_2.mol_list = None, None + with self.assertRaises(ReactionError): + ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), h2o_no_mol_2]) + + # 4. A species with neither coordinates nor a 2D graph: the check is skipped, as before. + empty = ARCSpecies(label='H2O', xyz=h2o_xyz) + empty.mol, empty.mol_list, empty.final_xyz, empty.initial_xyz, empty.conformers = None, None, None, None, list() + empty.most_stable_conformer, empty.cheap_conformer = None, None + self.assertEqual(empty.get_xyz(generate=False), None) + rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), empty]) + self.assertTrue(rxn.check_atom_balance()) + + # 5. A species with resonance structures: element counts are invariant under resonance. + c4h7 = ARCSpecies(label='C4H7', smiles='[CH2]C=CC') + self.assertGreater(len(c4h7.mol_list), 1) + for mol in c4h7.mol_list: + self.assertEqual(len(mol.atoms), len(c4h7.mol.atoms)) + rxn = ARCReaction(r_species=[c4h7, ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='C4H6', smiles='C=CC=C'), + ARCSpecies(label='H2O', smiles='O')]) + self.assertTrue(rxn.check_atom_balance()) + + # 6. A TS species, both with and without coordinates. + rxn = 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')]) + ts_xyz = {'symbols': ('C', 'H', 'H', 'H', 'H', 'O', 'H'), 'isotopes': (12, 1, 1, 1, 1, 16, 1), + 'coords': ((0.0, 0.0, 0.0), (1.09, 0.0, 0.0), (-0.36, 1.03, 0.0), (-0.36, -0.51, 0.89), + (-0.36, -0.51, -0.89), (0.0, 0.0, 2.5), (0.0, 0.0, 3.5))} + rxn.ts_species = ARCSpecies(label='TS', is_ts=True, xyz=ts_xyz) + self.assertTrue(rxn.check_atom_balance()) + rxn.ts_species = ARCSpecies(label='TS', is_ts=True) + self.assertTrue(rxn.check_atom_balance()) + + # 7. A TS species in a well: it reports no coordinates, and its graph is not used as a + # fallback, so the well is empty and the check is skipped even though the graph would + # have made the reaction unbalanced. This mirrors ARCSpecies.get_xyz(generate=True). + ts_in_well = ARCSpecies(label='ts_in_well', is_ts=True, smiles='CC') + self.assertIsNone(ts_in_well.get_xyz(generate=False)) + rxn = ARCReaction(r_species=[ts_in_well, ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='H2O', smiles='O')]) + self.assertTrue(rxn.check_atom_balance()) + + # 8. Species carrying coordinates use them, not the graph. + rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz), + ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz)], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz), + ARCSpecies(label='H2O', xyz=h2o_xyz)]) + self.assertTrue(rxn.check_atom_balance()) + + def test_get_atom_balance_entry(self): + """Test the _get_atom_balance_entry() function""" + # From a 2D graph, without generating a conformer. + ch4 = ARCSpecies(label='CH4', smiles='C') + entry = _get_atom_balance_entry(species=ch4) + self.assertEqual(sorted(line.split()[0] for line in entry.splitlines()), ['C', 'H', 'H', 'H', 'H']) + self.assertIsNone(ch4.cheap_conformer) + + # From available coordinates. + ch4_with_xyz = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) + self.assertEqual(_get_atom_balance_entry(species=ch4_with_xyz), + xyz_to_str(ch4_with_xyz.get_xyz(generate=False))) + + # From ``mol_list`` when ``mol`` was cleared, mirroring the + # ``self.mol is not None or self.mol_list is not None`` guard of ARCSpecies.get_xyz(). + mol_list_only = ARCSpecies(label='CH4', smiles='C') + mol_list_only.mol = None + self.assertIsNotNone(mol_list_only.mol_list) + self.assertEqual(sorted(line.split()[0] for line in + _get_atom_balance_entry(species=mol_list_only).splitlines()), + ['C', 'H', 'H', 'H', 'H']) + # The same state on a monoatomic, which reaches check_atom_balance() through a whole + # reaction. A monoatomic is served from its coordinates rather than from its graph, since + # ARCSpecies populates final_xyz for one on construction, so only the element is asserted. + h_mol_list_only = ARCSpecies(label='H', smiles='[H]') + h_mol_list_only.mol = None + self.assertEqual([line.split()[0] for line in + _get_atom_balance_entry(species=h_mol_list_only).splitlines()], ['H']) + rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), h_mol_list_only], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2', smiles='[H][H]')]) + self.assertTrue(rxn.check_atom_balance()) + + # Neither coordinates nor a graph. + empty = ARCSpecies(label='CH4', smiles='C') + empty.mol, empty.mol_list = None, None + self.assertEqual(_get_atom_balance_entry(species=empty), '') + + # A TS without coordinates never falls back to its graph. + ts = ARCSpecies(label='ts_no_xyz', is_ts=True, smiles='C') + self.assertEqual(_get_atom_balance_entry(species=ts), '') + + def test_check_atom_balance_does_not_seed_coordinates(self): + """Test that constructing a reaction leaves its graph-only species without coordinates. + + check_atom_balance() used to call get_xyz(generate=True) on every species, so every + reaction species silently acquired a force field geometry on construction. The scheduler + reads exactly that state through get_xyz(generate=False) when deciding whether to run + conformer jobs or to go straight to freq/sp, so the absence of the side effect is a + contract rather than an implementation detail. + """ + r_species = [ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='OH', smiles='[OH]')] + p_species = [ARCSpecies(label='CH3', smiles='[CH3]'), ARCSpecies(label='H2O', smiles='O')] + rxn = ARCReaction(r_species=r_species, p_species=p_species) + self.assertTrue(rxn.check_atom_balance()) + for species in r_species + p_species: + self.assertIsNone(species.get_xyz(generate=False), msg=f'{species.label} got coordinates') + self.assertIsNone(species.cheap_conformer, msg=f'{species.label} got a cheap conformer') + self.assertEqual(species.conformers, list(), msg=f'{species.label} got conformers') + + # A monoatomic is the documented exception: ARCSpecies populates its final_xyz on + # construction, before any reaction exists, since an atom has nothing to optimize. + h = ARCSpecies(label='H', smiles='[H]') + self.assertIsNotNone(h.get_xyz(generate=False)) + def test_get_species_count(self): """Test the get_species_count() method""" rxn1 = ARCReaction(reactants=['CH4', 'OH', 'H2O'], products=['CH3', 'H2O', 'H2O']) diff --git a/arc/scheduler.py b/arc/scheduler.py index f606976396..520465c2ec 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -1348,8 +1348,15 @@ def run_conformer_jobs(self, labels: list[str] | None = None): log_info_printed = True if self.species_dict[label].force_field == 'cheap': # Just embed in RDKit and use MMFF94s for opt and energies. - if self.species_dict[label].initial_xyz is None: - self.species_dict[label].initial_xyz = self.species_dict[label].get_xyz() + # The geometry is handed to process_conformers() as the single conformer rather + # than assigned to initial_xyz here: that method only spawns jobs for a species + # whose geometry is still unknown, so assigning initial_xyz first silences it + # and a species reaching this branch with opt turned off never gets its + # freq/sp/rotor jobs. + if self.species_dict[label].initial_xyz is None and not self.species_dict[label].conformers: + xyz = self.species_dict[label].get_xyz() + if xyz is not None: + self.species_dict[label].conformers = [xyz] else: # Run the combinatorial method w/o fitting a force field. n_confs = self.n_confs if self.species_dict[label].multi_species is None else 1 diff --git a/arc/scheduler_test.py b/arc/scheduler_test.py index 7803e010c2..a0c22be4fe 100644 --- a/arc/scheduler_test.py +++ b/arc/scheduler_test.py @@ -175,6 +175,39 @@ def setUpClass(cls): job_types=cls.job_types2, ) + def test_run_conformer_jobs_cheap_force_field_hands_over_a_conformer(self): + """Test that the 'cheap' force field path leaves its geometry where process_conformers() looks for it. + + A graph-only species reaches this path once check_atom_balance() stopped seeding coordinates + onto reaction species. process_conformers() only spawns jobs while the geometry is still + unknown, so the geometry must arrive as a conformer rather than as initial_xyz; assigning + initial_xyz here instead silences that method and a species whose 'opt' job type is off + never gets its freq/sp/rotor jobs. The job spawning itself is not asserted here because it + requires testing=False, which submits jobs to a server. + """ + spc = ARCSpecies(label='propane_cheap', smiles='CCC') + spc.force_field = 'cheap' + self.assertIsNone(spc.get_xyz(generate=False)) + sched = Scheduler(project='project_test_cheap_ff', ess_settings=self.ess_settings, + species_list=[spc], composite_method=None, + conformer_opt_level=Level(repr=default_levels_of_theory['conformer']), + 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']), + scan_level=Level(repr=default_levels_of_theory['scan']), + ts_guess_level=Level(repr=default_levels_of_theory['ts_guesses']), + project_directory=os.path.join(ARC_PATH, 'Projects', 'project_test_cheap_ff'), + testing=True, job_types={'conf_opt': False, 'opt': False, 'freq': True, 'sp': True, + 'rotors': False, 'irc': False, 'fine': False}, + orbitals_level=default_levels_of_theory['orbitals'], adaptive_levels=None, + ) + sched.run_conformer_jobs(labels=['propane_cheap']) + self.assertEqual(len(sched.species_dict['propane_cheap'].conformers), 1) + self.assertIsNone(sched.species_dict['propane_cheap'].initial_xyz) + self.assertEqual(sorted(sched.species_dict['propane_cheap'].conformers[0]['symbols']), + sorted(('C', 'C', 'C') + ('H',) * 8)) + shutil.rmtree(os.path.join(ARC_PATH, 'Projects', 'project_test_cheap_ff'), ignore_errors=True) + def test_conformers(self): """Test the parse_conformer_energy() and determine_most_stable_conformer() methods""" label = 'methylamine' diff --git a/arc/species/species.py b/arc/species/species.py index 7bc636ce25..25f0f2eaa8 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -1223,22 +1223,32 @@ def get_cheap_conformer(self): """ Cheaply (limiting the number of possible conformers) get a reasonable conformer, this could very well not be the best (lowest energy) one. + + The 2D graph is taken from ``mol_list`` when ``mol`` is unset, which is the same fallback + that ``is_monoatomic()``, ``is_diatomic()`` and the ``generate`` branch of ``get_xyz()`` + already make. Without it a polyatomic species holding only a ``mol_list`` raises an + ``AttributeError`` here rather than getting a conformer. """ + mol = self.mol if self.mol is not None else (self.mol_list[0] if self.mol_list else None) + if mol is None: + logger.warning(f'Could not generate a cheap conformer for {self.label}, it has no 2D graph.') + self.cheap_conformer = None + return if self.is_monoatomic(): self.cheap_conformer = \ - conformers.generate_monoatomic_conformer(symbol=self.mol_list[0].atoms[0].element.symbol)['xyz'] + conformers.generate_monoatomic_conformer(symbol=mol.atoms[0].element.symbol)['xyz'] self.initial_xyz = self.final_xyz = self.cheap_conformer elif self.is_diatomic(): self.cheap_conformer = \ - conformers.generate_diatomic_conformer(symbol_1=self.mol_list[0].atoms[0].element.symbol, - symbol_2=self.mol_list[0].atoms[1].element.symbol, + conformers.generate_diatomic_conformer(symbol_1=mol.atoms[0].element.symbol, + symbol_2=mol.atoms[1].element.symbol, multiplicity=self.multiplicity)['xyz'] else: - num_confs = min(500, max(50, len(self.mol.atoms) * 3)) - rd_mol = conformers.embed_rdkit(label=self.label, mol=self.mol, num_confs=num_confs) + num_confs = min(500, max(50, len(mol.atoms) * 3)) + rd_mol = conformers.embed_rdkit(label=self.label, mol=mol, num_confs=num_confs) xyzs, energies = conformers.rdkit_force_field(label=self.label, rd_mol=rd_mol, - mol=self.mol, + mol=mol, num_confs=num_confs, force_field='MMFF94s', ) diff --git a/arc/species/species_test.py b/arc/species/species_test.py index 1abac6ca85..fa2db0a08e 100644 --- a/arc/species/species_test.py +++ b/arc/species/species_test.py @@ -831,6 +831,25 @@ def test_get_number_of_electrons_does_not_generate_a_conformer(self): self.assertIsNone(spc.cheap_conformer) self.assertEqual(spc.conformers, list()) + def test_get_cheap_conformer_from_mol_list_only(self): + """Test that a polyatomic species holding only a ``mol_list`` still gets a cheap conformer. + + ``is_monoatomic()``, ``is_diatomic()`` and the ``generate`` branch of ``get_xyz()`` all + accept this state, so the polyatomic branch must accept it too rather than raise. + """ + spc = ARCSpecies(label='propane', smiles='CCC', compute_thermo=False) + spc.mol_list = spc.mol_list or [spc.mol] + spc.mol = None + xyz = spc.get_xyz(generate=True) + self.assertIsNotNone(xyz) + self.assertEqual(sorted(xyz['symbols']), sorted(('C', 'C', 'C') + ('H',) * 8)) + + # Neither a graph nor coordinates: a warning, not an exception. + no_structure = ARCSpecies(label='propane', smiles='CCC', compute_thermo=False) + no_structure.mol, no_structure.mol_list = None, None + no_structure.get_cheap_conformer() + self.assertIsNone(no_structure.cheap_conformer) + def test_as_dict(self): """Test Species.as_dict()""" spc_dict = self.spc3.as_dict()