diff --git a/arc/job/adapters/ts/heuristics.py b/arc/job/adapters/ts/heuristics.py index e0d227e56e..2954c0dc01 100644 --- a/arc/job/adapters/ts/heuristics.py +++ b/arc/job/adapters/ts/heuristics.py @@ -830,26 +830,33 @@ def find_distant_neighbor(mol: 'Molecule', def are_h_abs_wells_reversed(rxn: ARCReaction, product_dict: dict, - ) -> tuple[bool, bool]: + ) -> tuple[bool, bool, bool]: """ Determine whether the reactants or the products in an H_Abstraction reaction are reversed relative to the RMG template: R(*1)-H(*2) + R(*3)j <=> R(*1)j + R(*3)-H(*2) + ``reactants_reversed`` is True when the atom labeled *2 belongs to the reaction's second + reactant. ``products_reversed`` is True when the product bearing the atom labeled *2 is the + reaction's first product. ``dict_products_reversed`` is True when the atom labeled *2 belongs + to the first molecule of ``product_dict['products']``. ``product_dict['r_label_map']`` holds + zero-based global indices into the reaction's concatenated reactants, and + ``product_dict['p_label_map']`` holds zero-based global indices into the concatenated + ``product_dict['products']``. Args: rxn (ARCReaction): The ARCReaction object. product_dict (dict): The product dictionary. Returns: - tuple[bool, bool]: reactants_reversed, products_reversed. + tuple[bool, bool, bool]: reactants_reversed, products_reversed, dict_products_reversed. """ r_star_2 = product_dict['r_label_map']['*2'] p_star_2 = product_dict['p_label_map']['*2'] - r_species, p_species = rxn.get_reactants_and_products(return_copies=True) - reactants_reversed = len(r_species[0].mol.atoms) < r_star_2 - products_reversed = len(product_dict['products'][0].atoms) >= p_star_2 + r_species, p_species = rxn.get_reactants_and_products(return_copies=False) + reactants_reversed = r_star_2 >= len(r_species[0].mol.atoms) + dict_products_reversed = p_star_2 < len(product_dict['products'][0].atoms) same_order_between_rxn_prods_and_dict_prods = p_species[0].is_isomorphic(product_dict['products'][0]) - products_reversed = products_reversed == same_order_between_rxn_prods_and_dict_prods - return reactants_reversed, products_reversed + products_reversed = dict_products_reversed == same_order_between_rxn_prods_and_dict_prods + return reactants_reversed, products_reversed, dict_products_reversed def h_abstraction(reaction: ARCReaction, @@ -860,6 +867,10 @@ def h_abstraction(reaction: ARCReaction, ) -> list[dict]: """ Generate TS guesses for reactions of the RMG ``H_Abstraction`` family. + Only the ``H_Abstraction`` entries of ``reaction.product_dicts`` are considered, and the + orientation of the reactant and product wells relative to the RMG template is resolved + separately for each of them. An entry whose *2 label does not resolve to a hydrogen atom of the + respective reactant is skipped. Args: reaction: An ARCReaction instance. @@ -871,20 +882,26 @@ def h_abstraction(reaction: ARCReaction, dihedral_increment (int, optional): The dihedral increment to use for B-H-A-C and D-B-H-C dihedral scans. Returns: list[dict] - Entries are Cartesian coordinates of TS guesses for all reactions. Returns an empty list if - ``reaction.product_dicts`` is empty (i.e., no matching H_Abstraction family products were identified). + Entries are Cartesian coordinates of TS guesses for all reactions. + An empty list is returned if none of ``reaction.product_dicts`` belongs to + the ``H_Abstraction`` family. """ xyz_guesses = list() dihedral_increment = dihedral_increment or DIHEDRAL_INCREMENT - if not reaction.product_dicts: + product_dicts = [product_dict for product_dict in reaction.product_dicts + if product_dict.get('family') == 'H_Abstraction'] + if not product_dicts: logger.warning(f'Could not generate H_Abstraction TS guesses for reaction {reaction}: ' - f'no product_dicts were identified for this reaction. Other TS search methods will be attempted.') + f'no H_Abstraction product dictionaries were identified for it ' + f'(it has {len(reaction.product_dicts)} product dictionaries). ' + f'Other TS search methods will be attempted.') 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: + for product_dict in product_dicts: # 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. + reactants_reversed, products_reversed, dict_prods_reversed = are_h_abs_wells_reversed( + rxn=reaction, product_dict=product_dict) reactants, products = reaction.get_reactants_and_products(return_copies=False) reactant = reactants[int(reactants_reversed)] # Get R(*1)-H(*2). reactant_2 = reactants[int(not reactants_reversed)] # Get R(*3)j. @@ -894,10 +911,13 @@ def h_abstraction(reaction: ARCReaction, # Don't modify dihedrals for an attacking H (or other linear radical) at a linear angle, C ~ A -- H1 - H2 -- H. dihedral_increment = 360 h1 = product_dict['r_label_map']['*2'] - if reactants_reversed and h1 >= len(reactants[0].mol.atoms): + if reactants_reversed: h1 -= len(reactants[0].mol.atoms) + if h1 >= len(r_mol.atoms) or not r_mol.atoms[h1].is_hydrogen(): + logger.warning(f'Not using a product dictionary of reaction {reaction} in which the atom labeled *2 ' + f'is not a hydrogen atom of {reactant.label}.') + continue h2 = product_dict['p_label_map']['*2'] - dict_prods_reversed = h2 < len(product_dict['products'][0].atoms) dict_product = product_dict['products'][int(not dict_prods_reversed)] # Get R(*3)-H(*2) from the product_dict. product_atom_map = map_two_species(spc_1=dict_product, spc_2=product) if h2 >= len(product_atom_map): diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py index ea7de2d8b5..a3559306ef 100644 --- a/arc/job/adapters/ts/heuristics_test.py +++ b/arc/job/adapters/ts/heuristics_test.py @@ -9,6 +9,7 @@ import itertools import os import shutil +import tempfile import unittest from arc.common import ARC_TESTING_PATH, almost_equal_coords @@ -22,6 +23,7 @@ get_modified_params_from_zmat_2, get_new_map_based_on_zmat_1, get_new_zmat_2_map, + h_abstraction, stretch_zmat_bond, get_main_reactant_and_water_from_hydrolysis_reaction, setup_zmat_indices, @@ -32,9 +34,11 @@ check_ts_bonds, h_abstraction, ) +from arc.molecule.molecule import Molecule 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.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 @@ -559,6 +563,56 @@ def test_h_abstraction_with_empty_product_dicts(self): rxn1.product_dicts = [] xyz_guesses = h_abstraction(reaction=rxn1) self.assertEqual(xyz_guesses, []) + def test_heuristics_for_h_abstraction_symmetric_reactants(self): + """OH + OH <=> H2O + O, an H-abstraction between two identical reactants. + + The two product dicts label the transferring hydrogen (*2) in different reactants, so each of + them resolves ``reactants_reversed`` differently, and ``h1`` indexes a hydrogen of the + reactant that carries it in both. + """ + o_triplet = ARCSpecies(label='O', smiles='[O]', xyz='O 0.0 0.0 0.0') + rxn = ARCReaction(r_species=[self.oh, self.oh], p_species=[self.h2o, o_triplet]) + self.assertEqual(rxn.family, 'H_Abstraction') + reactants = rxn.get_reactants_and_products(return_copies=False)[0] + n_atoms_r_0 = len(reactants[0].mol.atoms) + reactants_reversed_values = list() + for product_dict in rxn.product_dicts: + reactants_reversed = are_h_abs_wells_reversed(rxn, product_dict=product_dict)[0] + reactants_reversed_values.append(reactants_reversed) + h1 = product_dict['r_label_map']['*2'] + if reactants_reversed: + self.assertGreaterEqual(h1, n_atoms_r_0) + h1 -= n_atoms_r_0 + else: + self.assertLess(h1, n_atoms_r_0) + reactant = reactants[int(reactants_reversed)] + self.assertLess(h1, len(reactant.mol.atoms)) + self.assertTrue(reactant.mol.atoms[h1].is_hydrogen()) + self.assertEqual(sorted(reactants_reversed_values), [False, True]) + project_directory = tempfile.mkdtemp() + self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True) + heuristics = HeuristicsAdapter(job_type='tsg', + reactions=[rxn], + testing=True, + project='test', + project_directory=project_directory, + dihedral_increment=30, + ) + heuristics.execute_incore() + self.assertGreater(len(rxn.ts_species.ts_guesses), 0) + for ts_guess in rxn.ts_species.ts_guesses: + xyz = ts_guess.initial_xyz + o_indices = [i for i, symbol in enumerate(xyz['symbols']) if symbol == 'O'] + h_indices = [i for i, symbol in enumerate(xyz['symbols']) if symbol == 'H'] + self.assertEqual(len(o_indices), 2) + self.assertEqual(len(h_indices), 2) + transferred_h = min(h_indices, + key=lambda i: sum(calculate_distance(xyz, [o, i]) for o in o_indices)) + for o in o_indices: + self.assertGreater(calculate_distance(xyz, [o, transferred_h]), 1.0) + self.assertLess(calculate_distance(xyz, [o, transferred_h]), 1.4) + self.assertGreater(calculate_angle(xyz, [o_indices[0], transferred_h, o_indices[1]]), 170) + self.assertGreater(calculate_distance(xyz, o_indices), 2.2) def test_heuristics_for_h_abstraction_2(self): # C3H8 + HO2 <=> C3H7 + H2O2 @@ -1998,7 +2052,7 @@ def _check_h_abs_wells_reversed(self, rxn: 'ARCReaction', expected_r_reversed: b consider_arc_families=False, discover_own_reverse_rxns_in_reverse=False, ) - r_reversed, p_reversed = are_h_abs_wells_reversed(rxn, product_dict=product_dicts[0]) + r_reversed, p_reversed, _ = are_h_abs_wells_reversed(rxn, product_dict=product_dicts[0]) self.assertEqual(r_reversed, expected_r_reversed) self.assertEqual(p_reversed, expected_p_reversed) @@ -2050,6 +2104,237 @@ def test_are_h_abs_wells_reversed_butyrate(self): ARCSpecies(label='H2O', smiles='O')]) self._check_h_abs_wells_reversed(rxn, expected_r_reversed=False, expected_p_reversed=False) + def _get_h_abs_product_dicts(self, rxn: 'ARCReaction') -> list: + """Get the H_Abstraction family product dicts of ``rxn``.""" + return get_reaction_family_products(rxn=rxn, + rmg_family_set=[rxn.family], + consider_rmg_families=True, + consider_arc_families=False, + discover_own_reverse_rxns_in_reverse=False, + ) + + def test_are_h_abs_wells_reversed_star_2_is_the_first_atom_of_the_second_reactant(self): + """CH3 + H2 <=> CH4 + H — *2 is the first atom of the second reactant. + + One of the two product dicts labels the first H of H2 as *2, i.e. r_label_map['*2'] equals + the number of atoms in the first reactant. The transferring H is then in the second + reactant, so reactants_reversed must be True, and h_abstraction must index h1 into H2 + rather than into CH3. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2', smiles='[H][H]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='H', smiles='[H]')]) + self.assertEqual(rxn.family, 'H_Abstraction') + product_dicts = self._get_h_abs_product_dicts(rxn) + n_atoms_r_0 = len(rxn.r_species[0].mol.atoms) + boundary_dicts = [p_dict for p_dict in product_dicts if p_dict['r_label_map']['*2'] == n_atoms_r_0] + self.assertEqual(len(boundary_dicts), 1) + reactants = rxn.get_reactants_and_products(return_copies=False)[0] + for product_dict in product_dicts: + reactants_reversed, _, _ = are_h_abs_wells_reversed(rxn, product_dict=product_dict) + reactant = reactants[int(reactants_reversed)] + h1 = product_dict['r_label_map']['*2'] + if reactants_reversed: + h1 -= len(reactants[0].mol.atoms) + self.assertLess(h1, len(reactant.mol.atoms)) + self.assertTrue(reactant.mol.atoms[h1].is_hydrogen()) + self.assertTrue(are_h_abs_wells_reversed(rxn, product_dict=boundary_dicts[0])[0]) + rxn.product_dicts = product_dicts + xyz_guesses = h_abstraction(reaction=rxn, dihedral_increment=120) + self.assertGreater(len(xyz_guesses), 0) + + def test_are_h_abs_wells_reversed_star_2_is_the_first_atom_of_the_second_dict_product(self): + """C2H6 + OH <=> C2H5 + H2O with the transferred H first in the second dict product. + + p_label_map['*2'] is a global index into the concatenated dict products, so a value equal + to the number of atoms of the first dict product already belongs to the second one, making + the H bearing dict product the second one. The dict products are constructed here with the + transferred H of the water first. The reaction product paired with the H bearing dict + product must be the water in both product orders. + """ + h2o_h_first = Molecule().from_adjacency_list("""1 H u0 p0 c0 {2,S} +2 O u0 p2 c0 {1,S} {3,S} +3 H u0 p0 c0 {2,S} +""") + c2h5, h2o = ARCSpecies(label='C2H5', smiles='[CH2]C'), ARCSpecies(label='H2O', smiles='O') + dict_products = [c2h5.mol.copy(deep=True), h2o_h_first] + p_star_2 = len(dict_products[0].atoms) + self.assertTrue(dict_products[1].atoms[0].is_hydrogen()) + rxn = ARCReaction(r_species=[ARCSpecies(label='C2H6', smiles='CC'), + ARCSpecies(label='OH', smiles='[OH]')], + p_species=[c2h5, h2o]) + self.assertEqual(rxn.family, 'H_Abstraction') + product_dict = self._get_h_abs_product_dicts(rxn)[0] + product_dict['products'] = dict_products + product_dict['p_label_map'] = {'*2': p_star_2} + _, products_reversed, dict_products_reversed = are_h_abs_wells_reversed(rxn, product_dict=product_dict) + self.assertFalse(dict_products_reversed) + self.assertFalse(products_reversed) + self.assertEqual(rxn.p_species[int(not products_reversed)].mol.get_formula(), 'H2O') + rxn_reversed_products = ARCReaction(r_species=[ARCSpecies(label='C2H6', smiles='CC'), + ARCSpecies(label='OH', smiles='[OH]')], + p_species=[h2o, c2h5]) + product_dict = self._get_h_abs_product_dicts(rxn_reversed_products)[0] + product_dict['products'] = dict_products + product_dict['p_label_map'] = {'*2': p_star_2} + _, products_reversed, dict_products_reversed = are_h_abs_wells_reversed(rxn_reversed_products, + product_dict=product_dict) + self.assertFalse(dict_products_reversed) + self.assertTrue(products_reversed) + self.assertEqual(rxn_reversed_products.p_species[int(not products_reversed)].mol.get_formula(), 'H2O') + + def test_are_h_abs_wells_reversed_with_a_charge_separated_perception(self): + """isoxazol-5-yl + propyne <=> isoxazole + propargyl, with a charge separated re-perception. + + Perceiving isoxazole from its Cartesian coordinates may yield the charge separated aromatic + [n-]1ccc[o+]1 rather than the neutral c1ccno1, and copying an ARCSpecies that carries + coordinates re-perceives its 2D graph. Here that re-perception is forced deterministically: + the orientation must still be resolved against the neutral graph the reaction holds, so that + the product paired with each dict product keeps its molecular formula. + """ + isoxazole_charge_separated_adjlist = """1 C u0 p0 c0 {2,S} {5,D} {6,S} +2 C u0 p0 c0 {1,S} {3,D} {7,S} +3 C u0 p0 c0 {2,D} {4,S} {8,S} +4 N u0 p2 c-1 {3,S} {5,S} +5 O u0 p1 c+1 {1,D} {4,S} +6 H u0 p0 c0 {1,S} +7 H u0 p0 c0 {2,S} +8 H u0 p0 c0 {3,S} +""" + isoxazolyl = ARCSpecies(label='isoxazolyl', smiles='O1[C]=CC=N1', multiplicity=2) + isoxazolyl.final_xyz = str_to_xyz("""O 1.62581768 -0.29842095 -0.47814520 +C 0.65307399 -1.11584355 -0.97117760 +C -0.57247729 -0.54346433 -0.73207366 +C -0.27021319 0.65775554 -0.07029399 +N 1.04374329 0.82411794 0.09216569 +H -1.53760759 -0.94291122 -1.00026412 +H -0.94233690 1.41876655 0.30089417""") + propyne = ARCSpecies(label='propyne', smiles='C#CC') + propyne.final_xyz = str_to_xyz("""C 1.69667781 -0.18348962 0.34176455 +C 0.50414063 -0.05452099 0.28779005 +C -0.94919940 0.10265252 0.22201144 +H 2.75513851 -0.29795839 0.38967067 +H -1.44299398 -0.87115669 0.15094730 +H -1.23736098 0.69413790 -0.65222854 +H -1.32640259 0.61033527 1.11485159""") + isoxazole = ARCSpecies(label='isoxazole', smiles='c1ccno1') + isoxazole.final_xyz = str_to_xyz("""C -1.09143635 -0.08868244 -0.00645500 +C -0.03176118 0.77687168 -0.08273051 +C 1.09337013 -0.05945157 0.02217020 +N 0.75345117 -1.34403600 0.15327038 +O -0.63329583 -1.35807858 0.13462360 +H -2.16536354 0.01990666 -0.03354656 +H -0.06648136 1.84914961 -0.19678701 +H 2.14151697 0.20432063 0.00945491""") + propargyl = ARCSpecies(label='propargyl', smiles='C#C[CH2]', multiplicity=2) + propargyl.final_xyz = str_to_xyz("""C 1.49185539 0.01644235 0.23786746 +C 0.30713117 0.00338503 0.04897028 +C -1.08917296 -0.01200422 -0.17366215 +H 2.54370249 0.02803524 0.40557840 +H -1.61674460 0.90953448 -0.38628957 +H -1.63677149 -0.94539289 -0.13246443""") + rxn = ARCReaction(r_species=[isoxazolyl, propyne], p_species=[isoxazole, propargyl]) + self.assertEqual(rxn.family, 'H_Abstraction') + self.assertEqual(rxn.p_species[0].mol.get_net_charge(), 0) + self.assertTrue(all(not atom.charge for atom in rxn.p_species[0].mol.atoms)) + product_dicts = self._get_h_abs_product_dicts(rxn) + original_mol_from_xyz = ARCSpecies.mol_from_xyz + + def charge_separated_mol_from_xyz(spc, xyz=None, get_cheap=False): + """Perceive isoxazole as its charge separated aromatic form, anything else as usual.""" + if spc.mol is not None and spc.mol.get_formula() == 'C3H3NO': + spc.mol = Molecule().from_adjacency_list(isoxazole_charge_separated_adjlist, + raise_atomtype_exception=False, + raise_charge_exception=False, + ) + return None + return original_mol_from_xyz(spc, xyz=xyz, get_cheap=get_cheap) + + self.addCleanup(setattr, ARCSpecies, 'mol_from_xyz', original_mol_from_xyz) + ARCSpecies.mol_from_xyz = charge_separated_mol_from_xyz + self.assertTrue(any(atom.charge for atom in isoxazole.copy().mol.atoms)) + for product_dict in product_dicts: + _, products_reversed, dict_products_reversed = are_h_abs_wells_reversed(rxn, product_dict=product_dict) + product = rxn.get_reactants_and_products(return_copies=False)[1][int(not products_reversed)] + dict_product = product_dict['products'][int(not dict_products_reversed)] + self.assertEqual(product.mol.get_formula(), dict_product.get_formula()) + + def test_h_abstraction_without_product_dicts(self): + """H2 + O <=> H + OH with no product dicts at all.""" + rxn = ARCReaction(r_species=[self.h2, self.o], p_species=[self.h, self.oh]) + self.assertEqual(rxn.family, 'H_Abstraction') + rxn.product_dicts = list() + with self.assertLogs('arc', level='WARNING') as log: + xyz_guesses = h_abstraction(reaction=rxn) + self.assertEqual(xyz_guesses, []) + warning = '\n'.join(log.output) + self.assertIn('H_Abstraction', warning) + self.assertIn(rxn.label, warning) + + def test_h_abstraction_with_product_dicts_of_another_family(self): + """CH3 + C2H5 <=> CH4 + C2H4, a disproportionation labeled as an H_Abstraction. + + The family attribute is set explicitly and is not derived from the product dicts, so the + product dicts of this reaction belong to another family and none of them may be read as an + H_Abstraction label map. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='C2H5', smiles='[CH2]C')], + p_species=[ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='C2H4', smiles='C=C')], + family='H_Abstraction') + self.assertEqual(rxn.family, 'H_Abstraction') + families = [product_dict['family'] for product_dict in rxn.product_dicts] + self.assertGreater(len(families), 0) + self.assertNotIn('H_Abstraction', families) + with self.assertLogs('arc', level='WARNING') as log: + xyz_guesses = h_abstraction(reaction=rxn) + self.assertEqual(xyz_guesses, []) + warning = '\n'.join(log.output) + self.assertIn('H_Abstraction', warning) + self.assertIn(rxn.label, warning) + + def test_h_abstraction_ignores_product_dicts_of_another_family(self): + """CH3 + H2 <=> CH4 + H with a product dict of another family added to the reaction.""" + rxn = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2', smiles='[H][H]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='H', smiles='[H]')]) + self.assertEqual(rxn.family, 'H_Abstraction') + h_abs_product_dicts = self._get_h_abs_product_dicts(rxn) + rxn.product_dicts = h_abs_product_dicts + xyz_guesses = h_abstraction(reaction=rxn, dihedral_increment=120) + self.assertGreater(len(xyz_guesses), 0) + rxn.product_dicts = h_abs_product_dicts + [{'family': 'Disproportionation', + 'products': list(), + 'r_label_map': {'*1': 0}, + 'p_label_map': {'*1': 0}}] + self.assertEqual(len(h_abstraction(reaction=rxn, dihedral_increment=120)), len(xyz_guesses)) + + def test_h_abstraction_skips_a_product_dict_whose_star_2_is_not_a_hydrogen(self): + """CH3 + H2 <=> CH4 + H with a product dict whose *2 label points at the carbon of CH3.""" + rxn = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2', smiles='[H][H]')], + p_species=[ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='H', smiles='[H]')]) + self.assertEqual(rxn.family, 'H_Abstraction') + product_dicts = self._get_h_abs_product_dicts(rxn) + self.assertEqual(len(product_dicts), 2) + rxn.product_dicts = product_dicts[:1] + xyz_guesses = h_abstraction(reaction=rxn, dihedral_increment=120) + self.assertGreater(len(xyz_guesses), 0) + corrupted_product_dict = dict(product_dicts[1]) + corrupted_product_dict['r_label_map'] = dict(corrupted_product_dict['r_label_map']) + corrupted_product_dict['r_label_map']['*2'] = 0 + self.assertFalse(rxn.r_species[0].mol.atoms[0].is_hydrogen()) + rxn.product_dicts = product_dicts[:1] + [corrupted_product_dict] + with self.assertLogs('arc', level='WARNING') as log: + self.assertEqual(len(h_abstraction(reaction=rxn, dihedral_increment=120)), len(xyz_guesses)) + warning = '\n'.join(log.output) + self.assertIn('*2', warning) + self.assertIn(rxn.label, warning) + def test_process_hydrolysis_reaction(self): """Test the process_hydrolysis_reaction() function.""" acetamide = self.acetamide diff --git a/arc/mapping/engine.py b/arc/mapping/engine.py index b67c90a0f4..d3f71794c6 100644 --- a/arc/mapping/engine.py +++ b/arc/mapping/engine.py @@ -21,8 +21,9 @@ from arc.species import ARCSpecies from arc.species.conformers import determine_chirality from arc.species.converter import compare_confs, sort_xyz_using_indices, xyz_from_data -from arc.species.vectors import (apply_rodrigues_rotation, calculate_dihedral_angle, get_angle, get_delta_angle, - get_perpendicular_axes, get_vector, get_vector_length, unit_vector) +from arc.species.vectors import (apply_rodrigues_rotation, calculate_dihedral_angle, + get_angle, get_delta_angle, get_perpendicular_axes, get_vector, + get_vector_length, is_torsion_linear, unit_vector) from arc.species.zmat import TOL_180 if TYPE_CHECKING: @@ -661,6 +662,9 @@ def get_backbone_dihedral_angles(spc_1: ARCSpecies, ) -> list[dict[str, float | list[int]]]: """ Determine the dihedral angles of the backbone torsions of two backbone mapped species. + Torsions that span a linear segment (e.g. a cumulene/ketene O=C=C backbone) are skipped: + their dihedral angle is geometrically undefined and ``ARCSpecies.set_dihedral()`` no-ops on + them, so feeding them into the backbone alignment loop only generates repeated log noise. The output has the following format:: torsions = [{'torsion 1': [0, 1, 2, 3], # The first torsion in terms of species 1's indices. @@ -689,8 +693,10 @@ def get_backbone_dihedral_angles(spc_1: ARCSpecies, if spc_1.mol.atoms[torsion_1[0]].is_non_hydrogen() \ and spc_1.mol.atoms[torsion_1[3]].is_non_hydrogen(): # This is not a "terminal" torsion. + torsion_2 = [backbone_map[t_1] for t_1 in torsion_1] + if is_torsion_linear(spc_1.get_xyz(), torsion_1) or is_torsion_linear(spc_2.get_xyz(), torsion_2): + continue for rotor_dict_2 in spc_2.rotors_dict.values(): - torsion_2 = [backbone_map[t_1] for t_1 in torsion_1] if all(pivot_2 in [torsion_2[1], torsion_2[2]] for pivot_2 in [rotor_dict_2['torsion'][1], rotor_dict_2['torsion'][2]]): torsions.append({'torsion 1': torsion_1, diff --git a/arc/mapping/engine_test.py b/arc/mapping/engine_test.py index 0c7e22c465..26944ef992 100644 --- a/arc/mapping/engine_test.py +++ b/arc/mapping/engine_test.py @@ -1266,6 +1266,33 @@ def test_get_backbone_dihedral_angles(self): self.assertAlmostEqual(torsions[0]['angle 1'], 67.81049913527622) self.assertAlmostEqual(torsions[0]['angle 2'], 174.65228274664804) + def test_is_torsion_linear(self): + """Test the is_torsion_linear() function. + + 2-pentyne (CC#CCC) has a collinear C1#C2-C3 alkyne segment, so the torsion [1, 2, 3, 4] + spans a linear segment (its [1, 2, 3] triplet is ~180 degrees), while a torsion around a + normal single bond (the terminal ethyl rotor) does not. + """ + spc = ARCSpecies(label='2-pentyne', smiles='CC#CCC') + spc.determine_rotors() + xyz = spc.get_xyz() + self.assertTrue(engine.is_torsion_linear(xyz, [1, 2, 3, 4])) + self.assertFalse(engine.is_torsion_linear(xyz, [2, 3, 4, 10])) + + def test_get_backbone_dihedral_angles_skips_linear_segment(self): + """Test that get_backbone_dihedral_angles() skips backbone torsions that span a linear segment. + + 2-pentyne's only heavy-atom-terminated backbone torsion, [1, 2, 3, 4], spans a linear alkyne + segment; it must be filtered out so it is never fed into the set_dihedral() alignment loop. + """ + spc_1 = ARCSpecies(label='2-pentyne-a', smiles='CC#CCC') + spc_2 = ARCSpecies(label='2-pentyne-b', smiles='CC#CCC') + spc_1.determine_rotors() + spc_2.determine_rotors() + backbone_map = {i: i for i in range(len(spc_1.mol.atoms))} + torsions = engine.get_backbone_dihedral_angles(spc_1, spc_2, backbone_map=backbone_map) + self.assertNotIn([1, 2, 3, 4], [torsion_dict['torsion 1'] for torsion_dict in torsions]) + def test_map_lists(self): """Test the map_lists function.""" self.assertEqual(engine.map_lists([], []), {}) diff --git a/arc/species/species.py b/arc/species/species.py index a326e21d8d..8dec131960 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -16,7 +16,6 @@ dfs, get_logger, get_single_bond_length, - is_angle_linear, is_multiplicity_parity_valid, is_xyz_linear, read_yaml_file, @@ -55,7 +54,7 @@ kabsch, ) from arc.species.perceive import perceive_molecule_from_xyz, is_mol_valid -from arc.species.vectors import calculate_angle, calculate_distance, calculate_dihedral_angle +from arc.species.vectors import calculate_distance, calculate_dihedral_angle, is_torsion_linear logger = get_logger() @@ -1472,9 +1471,8 @@ def set_dihedral(self, raise ValueError('Cannot set dihedral without xyz') if deg_increment is not None: deg_abs = calculate_dihedral_angle(coords=xyz, torsion=torsion) + deg_increment - if is_angle_linear(calculate_angle(coords=xyz, atoms=torsion[:3], index=0)) \ - or is_angle_linear(calculate_angle(coords=xyz, atoms=torsion[1:], index=0)): - logger.warning(f'Cannot change a dihedral that contains a linear segment. Got torsion:{torsion}, xyz:\n{xyz}') + if is_torsion_linear(xyz=xyz, torsion=torsion): + logger.debug(f'Cannot change a dihedral that contains a linear segment. Got torsion:{torsion}, xyz:\n{xyz}') return None mol = self.mol if mol is None: diff --git a/arc/species/species_test.py b/arc/species/species_test.py index ef48c7efc3..0305479379 100644 --- a/arc/species/species_test.py +++ b/arc/species/species_test.py @@ -383,6 +383,19 @@ def test_get_xyz(self): self.assertIsInstance(xyz, str) self.assertEqual(xyz, xyz_to_str(expected_xyz)) + def test_set_dihedral_linear_segment(self): + """Test that set_dihedral() no-ops (returns None) without a WARNING when the torsion spans a linear segment. + + Methylketene (CC=C=O) has a collinear C=C=O cumulene segment, so the torsion [0, 1, 2, 3] + contains a linear segment (the [1, 2, 3] triplet is ~180 degrees). + """ + spc = ARCSpecies(label='methylketene', smiles='CC=C=O') + xyz = spc.get_xyz() + with self.assertNoLogs(logger='arc', level='WARNING'): + result = spc.set_dihedral(scan=[0, 1, 2, 3], index=0, deg_abs=90.0, + count=False, chk_rotor_list=False, xyz=xyz) + self.assertIsNone(result) + def test_conformers(self): """Test conformer generation""" self.spc1.conformers = list() diff --git a/arc/species/vectors.py b/arc/species/vectors.py index 3cd4a587df..155bbc3c97 100644 --- a/arc/species/vectors.py +++ b/arc/species/vectors.py @@ -5,7 +5,7 @@ import math import numpy as np -from arc.common import logger +from arc.common import is_angle_linear, logger from arc.exceptions import VectorsError from arc.molecule.molecule import Molecule from arc.species._zmat_kernels import lib as _ck, available as _ck_available @@ -576,3 +576,22 @@ def get_delta_angle(a1: float, a1 %= 360 a2 %= 360 return min(abs(a1 - a2), abs(a1 + 360 - a2), abs(a1 - a2 - 360)) + + +def is_torsion_linear(xyz: dict, + torsion: list[int], + ) -> bool: + """ + Determine whether a torsion spans a linear segment: a ~180 degree angle over either of the + two atom triplets it is composed of, in which case its dihedral angle is geometrically + undefined and cannot be set or measured. + + Args: + xyz (dict): The 3D coordinates. + torsion (list[int]): The 0-indexed torsion atom indices. + + Returns: + bool: Whether the torsion contains a linear segment. + """ + return is_angle_linear(calculate_angle(coords=xyz, atoms=torsion[:3], index=0)) \ + or is_angle_linear(calculate_angle(coords=xyz, atoms=torsion[1:], index=0)) diff --git a/arc/species/vectors_test.py b/arc/species/vectors_test.py index 8efe80134a..cc81f66b49 100644 --- a/arc/species/vectors_test.py +++ b/arc/species/vectors_test.py @@ -495,6 +495,21 @@ def test_get_delta_angle(self): self.assertEqual(vectors.get_delta_angle(372, 359), 13.0) self.assertAlmostEqual(vectors.get_delta_angle(730, 10.1), 0.1) + def test_is_torsion_linear(self): + """Test identifying a torsion that spans a linear segment.""" + bent = {'symbols': ('H', 'O', 'O', 'H'), 'isotopes': (1, 16, 16, 1), + 'coords': ((0.0, 0.9, 0.4), (0.0, 0.0, 0.0), (1.5, 0.0, 0.0), (1.5, 0.9, -0.4))} + self.assertFalse(vectors.is_torsion_linear(xyz=bent, torsion=[0, 1, 2, 3])) + + ketene = {'symbols': ('H', 'C', 'C', 'O'), 'isotopes': (1, 12, 12, 16), + 'coords': ((0.0, 1.1, 0.0), (0.0, 0.0, 0.0), (0.0, -1.3, 0.0), (0.0, -2.5, 0.0))} + self.assertTrue(vectors.is_torsion_linear(xyz=ketene, torsion=[0, 1, 2, 3])) + + first_triplet_linear = {'symbols': ('H', 'C', 'C', 'H'), 'isotopes': (1, 12, 12, 1), + 'coords': ((0.0, 0.0, 0.0), (1.2, 0.0, 0.0), (2.4, 0.0, 0.0), + (3.0, 0.9, 0.0))} + self.assertTrue(vectors.is_torsion_linear(xyz=first_triplet_linear, torsion=[0, 1, 2, 3])) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))