diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 70c1b6576f..c6bcd3fe6f 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -19,7 +19,6 @@ ) from arc.imports import settings from arc.species.converter import check_xyz_dict, displace_xyz, xyz_to_dmat -from arc.mapping.engine import get_atom_indices_of_labeled_atoms_in_an_rmg_reaction from arc.statmech.factory import statmech_factory if TYPE_CHECKING: diff --git a/arc/common.py b/arc/common.py index dadd413bdb..d28af52beb 100644 --- a/arc/common.py +++ b/arc/common.py @@ -19,7 +19,6 @@ import warnings import yaml from collections import deque -from itertools import chain from typing import TYPE_CHECKING, Any, List, Optional, Tuple, Union import numpy as np @@ -27,7 +26,6 @@ import qcelemental as qcel from arkane.ess import ess_factory, GaussianLog, MolproLog, OrcaLog, QChemLog, TeraChemLog -import rmgpy from rmgpy.exceptions import AtomTypeError, ILPSolutionError, ResonanceError from rmgpy.molecule.atomtype import ATOMTYPES from rmgpy.molecule.element import get_element @@ -40,9 +38,7 @@ if TYPE_CHECKING: - from rmgpy.reaction import Reaction from rmgpy.species import Species - from arc.reaction import ARCReaction logger = logging.getLogger('arc') @@ -1563,50 +1559,6 @@ def calc_rmsd(x: Union[list, np.array], return float(rmsd) -def _check_r_n_p_symbols_between_rmg_and_arc_rxns(arc_reaction: 'ARCReaction', - rmg_reactions: List['Reaction'], - ) -> bool: - """ - A helper function to check that atom symbols are in the correct order between an ARC reaction - and its corresponding RMG reactions generated by the get_rmg_reactions_from_arc_reaction() function. - Used internally for testing. - - Args: - arc_reaction (ARCReaction): The ARCReaction object to inspect. - rmg_reactions (List['Reaction']): Entries are RMG Reaction objects to inspect. - Could contain either Species or Molecule object as reactants/products. - - Returns: - bool: Whether atom symbols are in the same respective order. - """ - result = True - num_rs, num_ps = len(arc_reaction.r_species), len(arc_reaction.p_species) - arc_r_symbols = [atom.element.symbol for atom in chain(*tuple(arc_reaction.r_species[i].mol.atoms for i in range(num_rs)))] - arc_p_symbols = [atom.element.symbol for atom in chain(*tuple(arc_reaction.p_species[i].mol.atoms for i in range(num_ps)))] - for rmg_reaction in rmg_reactions: - rmg_r_symbols = [atom.element.symbol - for atom in chain(*tuple(rmg_reaction.reactants[i].atoms - if isinstance(rmg_reaction.reactants[i], Molecule) - else rmg_reaction.reactants[i].molecule[0].atoms - for i in range(num_rs)))] - rmg_p_symbols = [atom.element.symbol - for atom in chain(*tuple(rmg_reaction.products[i].atoms - if isinstance(rmg_reaction.products[i], Molecule) - else rmg_reaction.products[i].molecule[0].atoms - for i in range(num_ps)))] - if any(symbol_1 != symbol_2 for symbol_1, symbol_2 in zip(arc_r_symbols, rmg_r_symbols)): - print('\nDifferent element order in reactants between ARC and RMG:') # Don't modify to logging. - print(arc_r_symbols) - print(rmg_r_symbols) - result = False - if any(symbol_1 != symbol_2 for symbol_1, symbol_2 in zip(arc_p_symbols, rmg_p_symbols)): - print('\nDifferent element order in products between ARC and RMG:') - print(arc_p_symbols) - print(rmg_p_symbols) - result = False - return result - - def safe_copy_file(source: str, destination: str, wait: int = 10, diff --git a/arc/common_test.py b/arc/common_test.py index 95d8737d94..9f665f6281 100644 --- a/arc/common_test.py +++ b/arc/common_test.py @@ -21,9 +21,7 @@ import arc.common as common from arc.exceptions import InputError, SettingsError from arc.imports import settings -from arc.mapping.engine import get_rmg_reactions_from_arc_reaction import arc.species.converter as converter -from arc.reaction import ARCReaction from arc.species.species import ARCSpecies @@ -1268,13 +1266,6 @@ def test_sort_atoms_in_descending_label_order(self): mol.atoms[0].label = "a" self.assertIsNone(common.sort_atoms_in_descending_label_order(mol=mol)) - def test_check_r_n_p_symbols_between_rmg_and_arc_rxns(self): - """Test the _check_r_n_p_symbols_between_rmg_and_arc_rxns() function""" - arc_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')]) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=arc_rxn) - self.assertTrue(common._check_r_n_p_symbols_between_rmg_and_arc_rxns(arc_rxn, rmg_reactions)) - def test_almost_equal_coords(self): """Test the almost_equal_coords() function""" with self.assertRaises(TypeError): diff --git a/arc/family/__init__.py b/arc/family/__init__.py index c316592387..24e9eca9db 100644 --- a/arc/family/__init__.py +++ b/arc/family/__init__.py @@ -1 +1,3 @@ import arc.family.family +from arc.family.family import ReactionFamily +from arc.family.family import get_reaction_family_products diff --git a/arc/family/family.py b/arc/family/family.py index 40b074bffc..ac027ae5da 100644 --- a/arc/family/family.py +++ b/arc/family/family.py @@ -11,7 +11,6 @@ from arc.common import clean_text, generate_resonance_structures, get_logger from arc.imports import settings -from arc.species.converter import check_isomorphism if TYPE_CHECKING: from arc.species import ARCSpecies @@ -519,7 +518,7 @@ def check_product_isomorphism(products: List['Molecule'], Returns: bool: Whether the products are isomorphic to the species. """ - prods_a = [generate_resonance_structures(mol) or [mol] for mol in products] + prods_a = [generate_resonance_structures(mol.copy(deep=True)) or [mol.copy(deep=True)] for mol in products] prods_b = [spc.mol_list or [spc.mol] for spc in p_species] if len(prods_a) == 1: prod_a = prods_a[0] @@ -564,6 +563,7 @@ def get_all_families(rmg_family_set: Union[List[str], str] = 'default', ) -> List[str]: """ Get all available RMG and ARC families. + If ``rmg_family_set`` is a list of family labels and does not contain family sets, it will be returned as is. Args: rmg_family_set (Union[List[str], str], optional): The RMG family set to use. @@ -574,9 +574,11 @@ def get_all_families(rmg_family_set: Union[List[str], str] = 'default', List[str]: A list of all available families. """ rmg_family_set = rmg_family_set or 'default' + family_sets = get_rmg_recommended_family_sets() + if isinstance(rmg_family_set, list) and all(fam not in family_sets for fam in rmg_family_set): + return rmg_family_set rmg_families, arc_families = list(), list() if consider_rmg_families: - family_sets = get_rmg_recommended_family_sets() if not isinstance(rmg_families, list) and rmg_family_set not in list(family_sets) + ['all']: raise ValueError(f'Invalid RMG family set: {rmg_family_set}') if rmg_family_set == 'all': diff --git a/arc/family/family_test.py b/arc/family/family_test.py index 2bec560809..ff25e4d3e7 100644 --- a/arc/family/family_test.py +++ b/arc/family/family_test.py @@ -605,6 +605,8 @@ def test_get_all_families(self): families = get_all_families(consider_rmg_families=False) self.assertIsInstance(families, list) self.assertIn('hydrolysis', families) + families = get_all_families(rmg_family_set=['H_Abstraction']) + self.assertEqual(families, ['H_Abstraction']) def test_get_rmg_recommended_family_sets(self): """Test getting RMG recommended family sets""" @@ -638,12 +640,12 @@ def test_load(self): self.assertFalse(fam_2.own_reverse) self.assertEqual(fam_2.reactants, [['Root']]) self.assertEqual(fam_2.product_num, 2) - self.assertEqual(fam_2.entries, {'Root': """1 *3 R!H u0 {2,S} {3,[S,D]} -2 *4 R!H u0 {1,S} {4,[S,D]} -3 *2 R!H u0 {1,[S,D]} {5,[D,T,B]} -4 *5 R!H u0 {2,[S,D]} {6,S} -5 *1 R!H u0 {3,[D,T,B]} -6 *6 H u0 {4,S}"""}) + self.assertEqual(fam_2.entries, {'Root': """1 *3 R!H u0 {2,S} {3,[S,D]} +2 *4 R!H u0 {1,S} {4,[S,D]} +3 *2 R!H u0 {1,[S,D]} {5,[D,T,B]} +4 *5 R!H u0 {2,[S,D]} {6,S} +5 *1 R!H u0 {3,[D,T,B]} +6 *6 [H,Li] u0 {4,S}"""}) self.assertEqual(fam_2.actions, [['CHANGE_BOND', '*1', -1, '*2'], ['BREAK_BOND', '*5', 1, '*6'], ['BREAK_BOND', '*3', 1, '*4'], @@ -960,6 +962,17 @@ def test_get_isomorphic_subgraph(self): ) self.assertEqual(isomorphic_subgraph, {0: '*3', 4: '*1', 7: '*2'}) + # def test_order_species_list(self): + # """Test the order_species_list() function""" + # spc1 = ARCSpecies(label='spc1', smiles='C') + # spc2 = ARCSpecies(label='spc2', smiles='CC') + # ordered_species_list = order_species_list(species_list=[spc2, spc1], reference_species=[spc1, spc2]) + # self.assertEqual(ordered_species_list, [spc1, spc2]) + # ordered_species_list = order_species_list(species_list=[spc2, spc1], reference_species=[spc2, spc1]) + # self.assertEqual(ordered_species_list, [spc2, spc1]) + # ordered_species_list = order_species_list(species_list=[spc2.mol, spc1], reference_species=[spc2, spc1.mol]) + # self.assertEqual(ordered_species_list, [spc2.mol, spc1]) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/arc/job/adapters/ts/heuristics.py b/arc/job/adapters/ts/heuristics.py index 76da11d2c9..635659db39 100644 --- a/arc/job/adapters/ts/heuristics.py +++ b/arc/job/adapters/ts/heuristics.py @@ -21,24 +21,21 @@ import numpy as np -from rmgpy.exceptions import ActionError from rmgpy.molecule.molecule import Molecule -from rmgpy.reaction import Reaction -from rmgpy.species import Species from arkane.statmech import is_linear from arc.common import almost_equal_coords, get_logger, is_angle_linear, key_by_val +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 -from arc.mapping.engine import map_arc_rmg_species, map_two_species +from arc.mapping.engine import map_two_species from arc.species.species import ARCSpecies, TSGuess, colliding_atoms from arc.species.zmat import get_parameter_from_atom_indices, remove_1st_atom, up_param if TYPE_CHECKING: - from rmgpy.data.kinetics.family import KineticsFamily from arc.level import Level from arc.reaction import ARCReaction @@ -254,20 +251,6 @@ def execute_incore(self): charge=rxn.charge, multiplicity=rxn.multiplicity, ) - rxn.arc_species_from_rmg_reaction() - reactants, products = rxn.get_reactants_and_products(arc=True, return_copies=True) - reactant_mol_combinations = list(itertools.product(*list(reactant.mol_list for reactant in reactants))) - product_mol_combinations = list(itertools.product(*list(product.mol_list for product in products))) - reaction_list = list() - for reactants in list(reactant_mol_combinations): - for products in list(product_mol_combinations): - rxns = react(reactants=list(reactants), - products=list(products), - family=rxn.family, - arc_reaction=rxn, - ) - if rxns is not None: - reaction_list.extend(rxns) xyzs = list() tsg = None @@ -281,10 +264,7 @@ def execute_incore(self): # r1_stretch_, r2_stretch_, a2_ = 1.2, 1.2, 170 # general guesses tsg = TSGuess(method='Heuristics') tsg.tic() - xyzs = h_abstraction(arc_reaction=rxn, - rmg_reactions=reaction_list, - dihedral_increment=self.dihedral_increment, - ) + xyzs = h_abstraction(arc_reaction=rxn, dihedral_increment=self.dihedral_increment) tsg.tok() for method_index, xyz in enumerate(xyzs): @@ -787,106 +767,6 @@ def update_new_map_based_on_zmat_2(new_map: dict, return new_map -def react(reactants: List[Union[Molecule, Species]], - products: List[Union[Molecule, Species]], - family: 'KineticsFamily', - arc_reaction: 'ARCReaction', - ) -> Optional[List[Reaction]]: - """ - React molecules to give the requested products via an RMG family, - resulting in a reaction with RMG's atom labels for the reactants and products. - - Args: - reactants (List['Molecule']): Entries are Molecule instances of the reaction reactants. - products (List['Molecule']): Entries are Molecule instances of the reaction products. - family (KineticsFamily): The RMG reaction family instance. - arc_reaction (ARCReaction): The corresponding ARCReaction object instance. - - Returns: - Optional[Reaction]: An RMG Reaction instance with atom-labeled reactants and products. - """ - # Assure Molecule object instances: - reactant_mols, product_mols = list(), list() - for reactant in reactants: - reactant_mols.append(reactant.copy(deep=True) if isinstance(reactant, Molecule) - else reactant.molecule[0].copy(deep=True)) - for product in products: - product_mols.append(product.copy(deep=True) if isinstance(product, Molecule) - else product.molecule[0].copy(deep=True)) - reactants_copy, products_copy = [r.copy(deep=True) for r in reactant_mols], [p.copy(deep=True) for p in product_mols] - - try: - reactions = family.generate_reactions(reactants=reactants_copy, - products=products_copy, - prod_resonance=False, - delete_labels=False, - relabel_atoms=False, - ) - except (ActionError, ValueError): - return None - for reaction in reactions: - try: - family.add_atom_labels_for_reaction(reaction=reaction, - output_with_resonance=False, - save_order=True, - ) - except (ActionError, ValueError): - continue - - for i in range(len(reactions)): - r_map, p_map = map_arc_rmg_species(arc_reaction=arc_reaction, rmg_reaction=reactions[i], concatenate=False) - ordered_rmg_reactants = [reactions[i].reactants[r_map[j]] for j in range(len(reactions[i].reactants))] - ordered_rmg_products = [reactions[i].products[p_map[j]] for j in range(len(reactions[i].products))] - reactions[i].reactants = ordered_rmg_reactants - reactions[i].products = ordered_rmg_products - - # Re-map all molecules since atoms in the RMG products were shuffled (products re-created from the reactants). - output, label_dicts = list(), list() - for reaction in reactions: - for index in range(len(reaction.reactants)): - # The RMG molecule will get a random 3D conformer, don't consider chirality when mapping. - atom_map = map_two_species(spc_1=reactant_mols[index], - spc_2=reaction.reactants[index], - consider_chirality=False, - ) - if atom_map is None: - break - new_atoms_list = list() - for i in range(len(reactant_mols[index].atoms)): - reactant_mols[index].atoms[i].id = reaction.reactants[index].molecule[0].atoms[atom_map[i]].id - reactant_mols[index].atoms[i].label = reaction.reactants[index].molecule[0].atoms[atom_map[i]].label - new_atoms_list.append(reactant_mols[index].atoms[i]) - reaction.reactants[index].molecule[0].atoms = new_atoms_list - else: - for index in range(len(reaction.products)): - # The RMG molecule will get a random 3D conformer, don't consider chirality when mapping. - atom_map = map_two_species(spc_1=product_mols[index], - spc_2=reaction.products[index], - consider_chirality=False, - ) - if atom_map is None: - break - new_atoms_list = list() - for i in range(len(product_mols[index].atoms)): - product_mols[index].atoms[i].id = reaction.products[index].molecule[0].atoms[atom_map[i]].id - product_mols[index].atoms[i].label = reaction.products[index].molecule[0].atoms[atom_map[i]].label - new_atoms_list.append(product_mols[index].atoms[i]) - reaction.products[index].molecule[0].atoms = new_atoms_list - else: - if reaction is not None: - # Check that the RMG reaction atom labels are unique. - label_dict = dict() - for reactant, product in zip(reaction.reactants, reaction.products): - for s, spc in enumerate([reactant, product]): - for i, atom in enumerate(spc.molecule[0].atoms): - if atom.label: - label_dict[f'{"P" if s else "R"}_{atom.label}'] = i - if label_dict not in label_dicts: - label_dicts.append(label_dict) - output.append(reaction) - return output - - def find_distant_neighbor(rmg_mol: 'Molecule', start: int, ) -> Optional[int]: @@ -918,8 +798,27 @@ def find_distant_neighbor(rmg_mol: 'Molecule', # Family-specific heuristics functions: -def h_abstraction(arc_reaction: 'ARCReaction', - rmg_reactions: List['Reaction'], +def are_h_abs_wells_reversed(rxn: 'ARCReaction', + product_dict: dict, + ) -> Tuple[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) + + Args: + rxn (ARCReaction): The ARCReaction object. + product_dict (dict): The product dictionary. + + Returns: + Tuple[bool, bool]: reactants_reversed, products_reversed. + """ + star_2 = product_dict['label_map']['*2'] + reactants_reversed = len(rxn.r_species[0].mol.atoms) <= star_2 + products_reversed = len(rxn.p_species[0].mol.atoms) <= rxn.atom_map[star_2] + return reactants_reversed, products_reversed + + +def h_abstraction(reaction: 'ARCReaction', r1_stretch: float = 1.2, r2_stretch: float = 1.2, a2: float = 180, @@ -929,10 +828,7 @@ def h_abstraction(arc_reaction: 'ARCReaction', Generate TS guesses for reactions of the RMG ``H_Abstraction`` family. Args: - arc_reaction: An ARCReaction instance. - rmg_reactions: Entries are RMGReaction instances. The reactants and products attributes should not contain - resonance structures as only the first molecule is considered -- pass several Reaction entries - instead. Atoms must be labeled according to the RMG reaction family. + 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 @@ -943,15 +839,41 @@ def h_abstraction(arc_reaction: 'ARCReaction', Returns: List[dict] Entries are Cartesian coordinates of TS guesses for all reactions. """ - if not len(rmg_reactions): - logger.warning(f'Cannot generate TS guesses for {arc_reaction} without an RMG Reaction object instance.') - return list() - xyz_guesses = list() dihedral_increment = dihedral_increment or DIHEDRAL_INCREMENT + product_dicts = get_reaction_family_products(rxn=reaction, + rmg_family_set=[reaction.family], + consider_rmg_families=True, + consider_arc_families=False, + discover_own_reverse_rxns_in_reverse=False, + ) # todo: loop through product_dict instead of rmg family, finish are_h_abs_wells_reversed test + expected_products = [{'discovered_in_reverse': False, + 'family': 'H_Abstraction', + 'group_labels': ('Xrad_H', 'Y_rad'), + 'label_map': {'*1': 0, '*2': 1, '*3': 3}, + 'own_reverse': True, + 'products': [Molecule(smiles="N"), Molecule(smiles="[NH]")]}, + {'discovered_in_reverse': False, + 'family': 'H_Abstraction', + 'group_labels': ('Xrad_H', 'Y_rad'), + 'label_map': {'*1': 0, '*2': 2, '*3': 3}, + 'own_reverse': True, + 'products': [Molecule(smiles="N"), Molecule(smiles="[NH]")]}, + {'discovered_in_reverse': False, + 'family': 'H_Abstraction', + 'group_labels': ('Y_rad', 'Xrad_H'), + 'label_map': {'*1': 3, '*2': 4, '*3': 0}, + 'own_reverse': True, + 'products': [Molecule(smiles="[NH]"), Molecule(smiles="N")]}, + {'discovered_in_reverse': False, + 'family': 'H_Abstraction', + 'group_labels': ('Y_rad', 'Xrad_H'), + 'label_map': {'*1': 3, '*2': 5, '*3': 0}, + 'own_reverse': True, + 'products': [Molecule(smiles="N"), Molecule(smiles="[NH]")]}] + # Identify R1H and R2H in the "R1H + R2 <=> R1 + R2H" or "R2 + R1H <=> R2H + R1" reaction - # using the first RMG reaction; all other RMG reactions and the ARC reaction should have the same order. # The expected RMG atom labels are: R(*1)-H(*2) + R(*3)j <=> R(*1)j + R(*3)-H(*2). reactants_reversed, products_reversed = False, False for atom in rmg_reactions[0].reactants[1].molecule[0].atoms: @@ -963,16 +885,21 @@ def h_abstraction(arc_reaction: 'ARCReaction', products_reversed = True break - arc_reactants, arc_products = arc_reaction.get_reactants_and_products(arc=True, return_copies=False) - arc_reactant = arc_reactants[int(reactants_reversed)] # Get R(*1)-H(*2). - arc_product = arc_products[int(not products_reversed)] # Get R(*3)-H(*2). + reactants, products = reaction.get_reactants_and_products(arc=True, return_copies=False) + reactant = reactants[int(reactants_reversed)] # Get R(*1)-H(*2). + product = products[int(not products_reversed)] # Get R(*3)-H(*2). - if any([is_linear(coordinates=np.array(arc_reactant.get_xyz()['coords'])), - is_linear(coordinates=np.array(arc_product.get_xyz()['coords']))]) and is_angle_linear(a2): + if any([is_linear(coordinates=np.array(reactant.get_xyz()['coords'])), + is_linear(coordinates=np.array(product.get_xyz()['coords']))]) and is_angle_linear(a2): # Don't modify dihedrals for an attacking H (or other linear radical) at a linear angle, C ~ A -- H1 - H2 -- H. dihedral_increment = 360 - for rmg_reaction in rmg_reactions: + for product_dict in product_dicts: + + reactants_reversed, products_reversed = are_h_abs_wells_reversed(rxn=reaction, product_dict=product_dict) + + + rmg_reactant_mol = rmg_reaction.reactants[int(reactants_reversed)].molecule[0] rmg_product_mol = rmg_reaction.products[int(not products_reversed)].molecule[0] h1 = rmg_reactant_mol.atoms.index([atom for atom in rmg_reactant_mol.atoms if atom.label == '*2'][0]) @@ -1002,11 +929,11 @@ def h_abstraction(arc_reaction: 'ARCReaction', xyz_guess = None try: xyz_guess = combine_coordinates_with_redundant_atoms( - xyz_1=arc_reactant.get_xyz(), - xyz_2=arc_product.get_xyz(), + xyz_1=reactant.get_xyz(), + xyz_2=product.get_xyz(), mol_1=rmg_reactant_mol, mol_2=rmg_product_mol, - reactant_2=arc_reaction.get_reactants_and_products()[0][int(not reactants_reversed)], + reactant_2=reaction.get_reactants_and_products()[0][int(not reactants_reversed)], h1=h1, h2=h2, c=c, diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py index 5d01d75650..93410cc8dd 100644 --- a/arc/job/adapters/ts/heuristics_test.py +++ b/arc/job/adapters/ts/heuristics_test.py @@ -14,8 +14,10 @@ from rmgpy.reaction import Reaction from rmgpy.species import Species -from arc.common import ARC_PATH, _check_r_n_p_symbols_between_rmg_and_arc_rxns, almost_equal_coords +from arc.common import ARC_PATH, almost_equal_coords +from arc.family import get_reaction_family_products from arc.job.adapters.ts.heuristics import (HeuristicsAdapter, + are_h_abs_wells_reversed, combine_coordinates_with_redundant_atoms, determine_glue_params, find_distant_neighbor, @@ -23,7 +25,6 @@ get_modified_params_from_zmat_2, get_new_map_based_on_zmat_1, get_new_zmat_2_map, - react, stretch_zmat_bond, ) from arc.reaction import ARCReaction @@ -1169,49 +1170,6 @@ def test_get_new_map_based_on_zmat_1(self): self.assertEqual(new_map, {0: 23, 1: 18, 2: 17, 3: 19, 4: 'X25', 5: 16, 6: 'X26', 7: 20, 8: 'X27', 9: 21, 10: 22, 11: 'X28', 12: 24}) # +16 - def test_react(self): - """Test the react() function and specifically that atom order in kept.""" - rxn_1 = ARCReaction(r_species=[ARCSpecies(label='C2H6', smiles='CC', xyz=self.c2h6_xyz), - ARCSpecies(label='CCOOj', smiles='CCO[O]', xyz=self.ccooj_xyz)], - p_species=[ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz), - ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz)]) - reactants, products = rxn_1.get_reactants_and_products(arc=True) - reactant_mol_combinations = list(itertools.product(*list(reactant.mol_list for reactant in reactants))) - product_mol_combinations = list(itertools.product(*list(product.mol_list for product in products))) - reactants = list(reactant_mol_combinations)[0] - products = list(product_mol_combinations)[0] - rmg_reactions = react(reactants=list(reactants), - products=list(products), - family=rxn_1.family, - arc_reaction=rxn_1, - ) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_1, rmg_reactions)) - - rxn_2 = ARCReaction(reactants=['H2N2(T)', 'N2H4'], products=['N2H3', 'N2H3'], - r_species=[ARCSpecies(label='H2N2(T)', smiles='[N]N'), - ARCSpecies(label='N2H4', smiles='NN')], - p_species=[ARCSpecies(label='N2H3', smiles='[NH]N')]) - self.assertEqual(rxn_2.family, 'H_Abstraction') - reactants, products = rxn_2.get_reactants_and_products(arc=False) - reactant_labels = [atom.label for atom in reactants[0].molecule[0].atoms if atom.label] \ - + [atom.label for atom in reactants[1].molecule[0].atoms if atom.label] - product_labels = [atom.label for atom in products[0].molecule[0].atoms if atom.label] \ - + [atom.label for atom in products[1].molecule[0].atoms if atom.label] - self.assertEqual(reactant_labels, list()) - self.assertEqual(product_labels, list()) - rmg_reactions = react(reactants=list(reactants), - products=list(products), - family=rxn_2.family, - arc_reaction=rxn_2, - ) - reactant_labels = [atom.label for atom in rmg_reactions[0].reactants[0].molecule[0].atoms if atom.label] \ - + [atom.label for atom in rmg_reactions[0].reactants[1].molecule[0].atoms if atom.label] - product_labels = [atom.label for atom in rmg_reactions[0].products[0].molecule[0].atoms if atom.label] \ - + [atom.label for atom in rmg_reactions[0].products[1].molecule[0].atoms if atom.label] - for label in ['*1', '*2', '*3']: - self.assertIn(label, reactant_labels) - self.assertIn(label, product_labels) - def test_generate_the_two_constrained_zmats(self): """Test the generate_the_two_constrained_zmats() function.""" zmat_1, zmat_2 = generate_the_two_constrained_zmats(xyz_1=self.ccooh_xyz, @@ -1701,6 +1659,29 @@ def test_find_distant_neighbor(self): self.assertEqual(find_distant_neighbor(rmg_mol=mol_3, start=2), 4) self.assertEqual(find_distant_neighbor(rmg_mol=mol_3, start=2), 4) + def test_are_h_abs_wells_reversed(self): + """Test the are_h_abs_wells_reversed() function.""" + rxn_1 = ARCReaction(r_species=[ARCSpecies(label='C2H6', smiles='CC'), ARCSpecies(label='OH', smiles='[OH]')], + p_species=[ARCSpecies(label='C2H5', smiles='[CH2]C'), ARCSpecies(label='H2O', smiles='O')]) + rxn_2 = ARCReaction(r_species=[ARCSpecies(label='OH', smiles='[OH]'), ARCSpecies(label='C2H6', smiles='CC')], # r reversed + p_species=[ARCSpecies(label='C2H5', smiles='[CH2]C'), ARCSpecies(label='H2O', smiles='O')]) + rxn_3 = ARCReaction(r_species=[ARCSpecies(label='C2H6', smiles='CC'), ARCSpecies(label='OH', smiles='[OH]')], # p reversed + p_species=[ARCSpecies(label='H2O', smiles='O'), ARCSpecies(label='C2H5', smiles='[CH2]C')]) + rxn_4 = ARCReaction(r_species=[ARCSpecies(label='OH', smiles='[OH]'), ARCSpecies(label='C2H6', smiles='CC')], # r and p reversed + p_species=[ARCSpecies(label='H2O', smiles='O'), ARCSpecies(label='C2H5', smiles='[CH2]C')]) + + product_dicts = get_reaction_family_products(rxn=rxn_1, + rmg_family_set=[rxn_1.family], + consider_rmg_families=True, + consider_arc_families=False, + discover_own_reverse_rxns_in_reverse=False, + ) + r_reversed, p_reversed = are_h_abs_wells_reversed(rxn_1, product_dict=product_dicts[0]) + self.assertFalse(r_reversed) + self.assertFalse(p_reversed) + + + @classmethod def tearDownClass(cls): """ diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index bbdd176183..5df15faf27 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -7,15 +7,12 @@ 3) If the reaction is supported by RMG, it is sent to the driver. Else, it is mapped with map_general_rxn. """ -from typing import TYPE_CHECKING, List, Optional +from typing import TYPE_CHECKING, Dict, List, Optional -from arc.mapping.engine import (assign_labels_to_products, - are_adj_elements_in_agreement, +from arc.mapping.engine import (are_adj_elements_in_agreement, create_qc_mol, flip_map, fingerprint, - get_atom_indices_of_labeled_atoms_in_an_rmg_reaction, - get_rmg_reactions_from_arc_reaction, glue_maps, label_species_atoms, make_bond_changes, @@ -23,22 +20,22 @@ iterative_dfs, map_two_species, pairing_reactants_and_products_for_mapping, copy_species_list_for_mapping, - find_all_bdes, + find_all_breaking_bonds, cut_species_based_on_atom_indices, update_xyz, RESERVED_FINGERPRINT_KEYS,) from arc.common import logger +from arc.species.converter import check_molecule_list_order from rmgpy.exceptions import ActionError, AtomTypeError if TYPE_CHECKING: - from rmgpy.data.rmg import RMGDatabase + from rmgpy.molecule.molecule import Molecule from arc.reaction import ARCReaction def map_reaction(rxn: 'ARCReaction', backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, flip = False ) -> Optional[List[int]]: """ @@ -47,7 +44,6 @@ def map_reaction(rxn: 'ARCReaction', Args: rxn (ARCReaction): An ARCReaction object instance. backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. - db (RMGDatabase, optional): The RMG database instance. Returns: Optional[List[int]]: @@ -76,7 +72,6 @@ def map_reaction(rxn: 'ARCReaction', def map_general_rxn(rxn: 'ARCReaction', backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, ) -> Optional[List[int]]: """ Map a general reaction (one that was not categorized into a reaction family by RMG). @@ -85,7 +80,6 @@ def map_general_rxn(rxn: 'ARCReaction', Args: rxn (ARCReaction): An ARCReaction object instance. backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. - db (RMGDatabase, optional): The RMG database instance. Returns: Optional[List[int]]: @@ -121,7 +115,6 @@ def map_isomerization_reaction(rxn: 'ARCReaction', Args: rxn (ARCReaction): An ARCReaction object instance. backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. - db (RMGDatabase, optional): The RMG database instance. Returns: Optional[List[int]]: @@ -209,7 +202,6 @@ def map_isomerization_reaction(rxn: 'ARCReaction', def map_rxn(rxn: 'ARCReaction', backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, ) -> Optional[List[int]]: """ A wrapper function for mapping reaction, uses databases for mapping with the correct reaction family parameters. @@ -224,7 +216,6 @@ def map_rxn(rxn: 'ARCReaction', Args: rxn (ARCReaction): An ARCReaction object instance that belongs to the RMG H_Abstraction reaction family. backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. - db (RMGDatabase, optional): The RMG database instance. Returns: Optional[List[int]]: @@ -232,7 +223,7 @@ def map_rxn(rxn: 'ARCReaction', corresponding entry values are running atom indices of the products. """ # step 1: - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn, backend=backend) + rmg_reactions = 5 if not rmg_reactions: return None @@ -243,11 +234,11 @@ def map_rxn(rxn: 'ARCReaction', # step 2: assign_labels_to_products(rxn, p_label_dict) - #step 3: + # step 3: reactants, products = copy_species_list_for_mapping(rxn.r_species), copy_species_list_for_mapping(rxn.p_species) label_species_atoms(reactants), label_species_atoms(products) - r_bdes, p_bdes = find_all_bdes(rxn, r_label_dict, True), find_all_bdes(rxn, p_label_dict, False) + r_bdes, p_bdes = find_all_breaking_bonds(rxn, True), find_all_breaking_bonds(rxn, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) @@ -259,7 +250,7 @@ def map_rxn(rxn: 'ARCReaction', r_cuts, p_cuts = update_xyz(r_cuts), update_xyz(p_cuts) - #step 4: + # step 4: pairs_of_reactant_and_products = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts) if len(p_cuts): logger.error(f"Could not find isomorphism for scissored species: {[cut.mol.smiles for cut in p_cuts]}") @@ -267,5 +258,45 @@ def map_rxn(rxn: 'ARCReaction', # step 5: maps = map_pairs(pairs_of_reactant_and_products) - #step 6: + # step 6: return glue_maps(maps, pairs_of_reactant_and_products) + + +def convert_label_dict(label_dict: Dict[str, int], + reference_mol_list: List['Molecule'], + mol_list: List['Molecule'], + ) -> Optional[Dict[str, int]]: + """ + Convert the label dictionary to the correct atom indices in the reaction and reference molecules + + Args: + label_dict (Dict[str, int]): A dictionary of atom labels (e.g., '*1') to atom indices. + reference_mol_list (List[Molecule]): The list of molecules to which label_dict values refer. + mol_list (List[Molecule]): The list of molecules to which label_dict values should be converted. + + Returns: + Dict[str, int]: The converted label dictionary. + """ + print(f'original label_dict: {label_dict}') + if len(reference_mol_list) != len(mol_list): + raise ValueError(f'The number of reference molecules ({len(reference_mol_list)}) ' + f'does not match the number of molecules ({len(mol_list)}).') + if len(reference_mol_list) == 1: + atom_map = map_two_species(reference_mol_list[0], mol_list[0]) + if atom_map is None: + print(f'Could not map {reference_mol_list[0].to_smiles()} to {mol_list[0].to_smiles()}') + return None + return {label: atom_map[index] for label, index in label_dict.items()} + elif len(reference_mol_list) == 2: + ordered = check_molecule_list_order(mols_1=reference_mol_list, mols_2=mol_list) + print(f'ordered: {ordered}') + atom_map_1 = map_two_species(reference_mol_list[0], mol_list[0]) if ordered else map_two_species(reference_mol_list[1], mol_list[0]) + atom_map_2 = map_two_species(reference_mol_list[1], mol_list[1]) if ordered else map_two_species(reference_mol_list[0], mol_list[1]) + if atom_map_1 is None or atom_map_2 is None: + print(f'Could not map {reference_mol_list[0].to_smiles()} to {mol_list[0].to_smiles()} ' + f'or {reference_mol_list[1].to_smiles()} to {mol_list[1].to_smiles()}') + return None + atom_map = atom_map_1 + [index + len(atom_map_1) for index in atom_map_2] if ordered else \ + atom_map_2 + [index + len(atom_map_2) for index in atom_map_1] + print(f'atom_map: {atom_map}') + return {label: atom_map[index] for label, index in label_dict.items()} diff --git a/arc/mapping/driver_test.py b/arc/mapping/driver_test.py index dba1f085c2..f35e7156b1 100644 --- a/arc/mapping/driver_test.py +++ b/arc/mapping/driver_test.py @@ -11,6 +11,7 @@ from rmgpy.reaction import Reaction from rmgpy.species import Species +from arc.family import get_reaction_family_products from arc.mapping.driver import * from arc.reaction import ARCReaction from arc.mapping.engine import check_atom_map @@ -1057,6 +1058,77 @@ def test_map_isomerization_reaction(self): self.assertIn(atom_map[7], [6, 7]) self.assertIn(atom_map[8], [7, 8]) + def test_convert_label_dict(self): + """Test the convert_label_dict() function.""" + rxn_1 = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C'), ARCSpecies(label='O2', smiles='[O][O]')], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]'), ARCSpecies(label='HO2', smiles='O[O]')]) + products = get_reaction_family_products(rxn_1) + + self.assertEqual(products[0]['p_label_map'], {'*1': 3, '*2': 2, '*3': 1}) + p_label_dict = convert_label_dict(label_dict=products[0]['p_label_map'], + reference_mol_list=products[0]['products'], + mol_list=[spc.mol for spc in rxn_1.p_species]) + print(p_label_dict) + self.assertEqual(p_label_dict, {'*1': 3, '*2': 2, '*3': 1}) # did it do anyhting? looks like it didnt because the order should be different + + + # expected_products = [{'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 1, '*3': 5}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 1}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 1, '*3': 6}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 0}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 2, '*3': 5}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 1}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 2, '*3': 6}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 0}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 3, '*3': 5}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 1}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 3, '*3': 6}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 0}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 4, '*3': 5}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 1}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}, + # {'discovered_in_reverse': False, + # 'family': 'H_Abstraction', + # 'group_labels': ('X_H', 'Y_rad'), + # 'r_label_map': {'*1': 0, '*2': 4, '*3': 6}, + # 'p_label_map': {'*1': 3, '*2': 2, '*3': 0}, + # 'own_reverse': True, + # 'products': [Molecule(smiles="[O]O"), Molecule(smiles="[CH3]")]}] + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/arc/mapping/engine.py b/arc/mapping/engine.py index 1d8e706728..3f32dbaf14 100644 --- a/arc/mapping/engine.py +++ b/arc/mapping/engine.py @@ -19,7 +19,7 @@ from arc.common import convert_list_index_0_to_1, extremum_list, generate_resonance_structures, logger, key_by_val from arc.exceptions import SpeciesError -from arc.family import ReactionFamily +from arc.family import ReactionFamily, get_reaction_family_products from arc.species import ARCSpecies from arc.species.conformers import determine_chirality from arc.species.converter import compare_confs, sort_xyz_using_indices, translate_xyz, xyz_from_data, xyz_to_str @@ -735,7 +735,7 @@ def map_hydrogens(spc_1: ARCSpecies, return atom_map -def check_atom_map(rxn) -> bool: +def check_atom_map(rxn: 'ARCReaction') -> bool: """ A helper function for testing a reaction atom map. Tests that element symbols are ordered correctly. @@ -834,7 +834,8 @@ def flip_map(atom_map: Optional[List[int]]) -> Optional[List[int]]: def make_bond_changes(rxn: 'ARCReaction', r_cuts: List[ARCSpecies], - r_label_dict: Dict): + r_label_dict: Dict, + ) -> None: """ Makes the bond change before matching the reactants and products @@ -843,16 +844,17 @@ def make_bond_changes(rxn: 'ARCReaction', r_cuts: the cut products r_label_dict: the dictionary object the find the relevant location. """ - for action in rxn.family.forward_recipe.actions: - if action[0].lower() == "CHANGE_BOND".lower(): - indicies = r_label_dict[action[1]],r_label_dict[action[3]] + family = ReactionFamily(label=rxn.family) + for action in family.actions: + if action[0].lower() == "change_bond": + indices = r_label_dict[action[1]],r_label_dict[action[3]] for r_cut in r_cuts: - if indicies[0] in [int(atom.label) for atom in r_cut.mol.atoms] and indicies[1] in [int(atom.label) for atom in r_cut.mol.atoms]: + if indices[0] in [int(atom.label) for atom in r_cut.mol.atoms] and indices[1] in [int(atom.label) for atom in r_cut.mol.atoms]: atom1, atom2 = 0, 0 for atom in r_cut.mol.atoms: - if int(atom.label) == indicies[0]: + if int(atom.label) == indices[0]: atom1 = atom - elif int(atom.label) == indicies[1]: + elif int(atom.label) == indices[1]: atom2 = atom if atom1.radical_electrons == 0 and atom2.radical_electrons == 0: # Both atoms do not have any radicals, but their bond is changing. There probably is resonance, so this will not affect the isomorphism check. return @@ -945,14 +947,13 @@ def map_pairs(pairs: List[Tuple[ARCSpecies, ARCSpecies]]) -> List[List[int]]: Returns: List[List[int]]: A list of the mapped species """ - maps = list() for pair in pairs: maps.append(map_two_species(pair[0], pair[1])) return maps -def label_species_atoms(species: List['ARCSpecies']): +def label_species_atoms(species: List['ARCSpecies']) -> None: """ Adds the labels to the ``.mol.atoms`` properties of the species object. @@ -966,9 +967,11 @@ def label_species_atoms(species: List['ARCSpecies']): index += 1 -def glue_maps(maps: List[List[int]], pairs_of_reactant_and_products: List[Tuple[ARCSpecies, ARCSpecies]]) -> List[int]: +def glue_maps(maps: List[List[int]], + pairs_of_reactant_and_products: List[Tuple[ARCSpecies, ARCSpecies]], + ) -> List[int]: """ - a function that joins together the maps from the parts of the reaction. + Join the maps from the parts of a bimolecular reaction. Args: maps (List[List[int]]): The list of all maps of the isomorphic cuts. @@ -998,7 +1001,6 @@ def determine_bdes_on_spc_based_on_atom_labels(spc: "ARCSpecies", bde: Tuple[int Returns: bool: Whether the bde is based on the atom labels. """ - bde = convert_list_index_0_to_1(bde, direction=-1) index1, index2 = bde[0], bde[1] new_bde, atoms = list(), list() for index, atom in enumerate(spc.mol.atoms): @@ -1007,12 +1009,10 @@ def determine_bdes_on_spc_based_on_atom_labels(spc: "ARCSpecies", bde: Tuple[int atoms.append(atom) if len(new_bde) == 2: break - if len(new_bde) == 2 and atoms[1] in atoms[0].bonds.keys(): spc.bdes = [tuple(new_bde)] return True - else: - return False + return False def cut_species_based_on_atom_indices(species: List["ARCSpecies"], @@ -1032,7 +1032,7 @@ def cut_species_based_on_atom_indices(species: List["ARCSpecies"], Optional[List["ARCSpecies"]]: The species list input after the scission. """ if ref_species is not None: - bdes = translate_indices_based_on_ref_species(species, ref_species, bdes) + bdes = translate_bdes_based_on_ref_species(species, ref_species, bdes) if not bdes: return species for bde in bdes: @@ -1058,21 +1058,61 @@ def cut_species_based_on_atom_indices(species: List["ARCSpecies"], return species +def translate_bdes_based_on_ref_species(species: List["ARCSpecies"], + ref_species: List["ARCSpecies"], + bdes: List[Tuple[int, int]], + ) -> Optional[List[Tuple[int, int]]]: + """ + Translate a list of BDE indices based on a reference species list. + The given BDE indices refer to `ref_species`, and they'll be translated to refer to `species`. + + Args: + species (List[ARCSpecies]): The species list for which the indices should be translated. + ref_species (List[ARCSpecies]): The reference species list. + bdes (List[Tuple[int, int]]): The BDE indices to be translated. + + Returns: + Optional[List[Tuple[int, int]]]: The translated BDE indices, or None if translation fails. + """ + if not bdes: + return list() + all_indices = set() + for bde in bdes: + all_indices.update(bde) + sorted_indices = sorted(all_indices) + translated_indices = translate_indices_based_on_ref_species(species=species, + ref_species=ref_species, + indices=sorted_indices) + if translated_indices is None: + return None + index_translation_map = dict(zip(sorted_indices, translated_indices)) + new_bdes = list() + for bde in bdes: + a, b = bde + translated_a = index_translation_map.get(a) + translated_b = index_translation_map.get(b) + if translated_a is not None and translated_b is not None: + new_bdes.append(tuple(sorted([translated_a, translated_b]))) + else: + return None + return new_bdes + + def translate_indices_based_on_ref_species(species: List["ARCSpecies"], ref_species: List["ARCSpecies"], - bdes: List[Tuple[int, int]], - ) -> Optional[List[Tuple[int, int]]]: + indices: List[int], + ) -> Optional[List[int]]: """ - A function for translating the atom indices based on a reference species list. - The given bde indices refer to ``ref_species``, and they'll be translated to refer to ``species``. + Translate a list of atom indices based on a reference species list. + The given indices refer to `ref_species`, and they'll be translated to refer to `species`. Args: species (List[ARCSpecies]): The species list for which the indices should be translated. ref_species (List[ARCSpecies]): The reference species list. - bdes (List[Tuple[int, int]]): The BDE indices to be translated. + indices (List[int]): The list of atom indices to be translated. Returns: - Optional[List[Tuple[int, int]]]: The translated BDE indices. + Optional[List[int]]: The translated atom indices if all translations are successful; otherwise, None. """ visited_ref_species = list() species_map = dict() # maps ref species j to species i @@ -1083,42 +1123,43 @@ def translate_indices_based_on_ref_species(species: List["ARCSpecies"], visited_ref_species.append(j) species_map[j] = i index_map[j] = map_two_species(ref_spc, spc) + print(f'Found a match between species {i} and ref species {j}, map: {index_map[j]}') break - new_bdes = list() ref_spcs_lengths = [ref_spc.number_of_atoms for ref_spc in ref_species] accum_sum_ref_spcs_lengths = [sum(ref_spcs_lengths[:i+1]) for i in range(len(ref_spcs_lengths))] spcs_lengths = [spc.number_of_atoms for spc in species] accum_sum_spcs_lengths = [sum(spcs_lengths[:i+1]) for i in range(len(spcs_lengths))] - for bde in bdes: - a, b = bde - translated_bde = list() - for n in [a, b]: - found = False - for j, ref_len in enumerate(accum_sum_ref_spcs_lengths): - if n < ref_len: - atom_map = index_map[j] - i = species_map[j] - if atom_map is None or i is None: - return None - increment = accum_sum_spcs_lengths[i - 1] if i > 0 else 0 - translated_atom = atom_map[n - accum_sum_ref_spcs_lengths[j]] + increment - translated_bde.append(translated_atom) - found = True - break - if not found: - return None - if len(translated_bde) == 2: - new_bdes.append(tuple(translated_bde)) + + def translate_single_index(index: int) -> Optional[int]: + for ref_idx, ref_len in enumerate(accum_sum_ref_spcs_lengths): + if index < ref_len: + atom_map = index_map.get(ref_idx) + species_i = species_map.get(ref_idx) + if atom_map is None or species_i is None: + return None + increment = accum_sum_spcs_lengths[species_i - 1] if species_i > 0 else 0 + ref_start = accum_sum_ref_spcs_lengths[ref_idx - 1] if ref_idx > 0 else 0 + translated_atom = atom_map[index - ref_start] + increment + return translated_atom + return None + + new_indices = list() + for idx in indices: + translated_idx = translate_single_index(idx) + if translated_idx is not None: + new_indices.append(translated_idx) else: return None - return new_bdes + return new_indices def copy_species_list_for_mapping(species: List["ARCSpecies"]) -> List["ARCSpecies"]: """ A helper function for copying the species list for mapping. Also keeps the atom indices when copying. + Args: species (List[ARCSpecies]): The species list to be copied. + Returns: List[ARCSpecies]: The copied species list. """ @@ -1130,7 +1171,6 @@ def copy_species_list_for_mapping(species: List["ARCSpecies"]) -> List["ARCSpeci def find_all_breaking_bonds(rxn: "ARCReaction", - label_dict: Dict[str, int], r_direction: bool, ) -> Optional[List[Tuple[int, int]]]: """ @@ -1139,16 +1179,28 @@ def find_all_breaking_bonds(rxn: "ARCReaction", Args: rxn (ARCReaction): The reaction in question. - label_dict (Dict[str, int]): Keys are atom labels (e.g., '*1'), values are atom indices (0-indexed). r_direction (bool): Whether to consider the reactants direction (``True``) or the products direction (``False``). Returns: List[Tuple[int, int]]: Entries are tuples of the form (atom_index1, atom_index2) for each broken bond (1-indexed), representing the atom indices to be cut. """ + print(f'*** rxn.family: {rxn.family}') # todo: product bdes are not being translated correctly family = ReactionFamily(label=rxn.family) + product_dicts = get_reaction_family_products(rxn=rxn, rmg_family_set=[rxn.family]) + print(f'*** product_dicts: {product_dicts}') + label_dict = product_dicts[0]['r_label_map'] if r_direction else product_dicts[0]['p_label_map'] + print(f'Label dict: {label_dict}') breaking_bonds = list() for action in family.actions: if action[0].lower() == ("break_bond" if r_direction else "form_bond"): - breaking_bonds.append((label_dict[action[1]], label_dict[action[3]])) + breaking_bonds.append(tuple(sorted((label_dict[action[1]], label_dict[action[3]])))) + print(f'appended {breaking_bonds[-1]}') + print(f'original Breaking bonds: {breaking_bonds}') + if not r_direction: + breaking_bonds = translate_bdes_based_on_ref_species( + species=rxn.get_reactants_and_products()[1], + ref_species=[ARCSpecies(label=f'S{i}', mol=mol) for i, mol in enumerate(product_dicts[0]['products'])], + bdes=breaking_bonds) + print(f'translated Breaking bonds: {breaking_bonds}') return breaking_bonds diff --git a/arc/mapping/engine_test.py b/arc/mapping/engine_test.py index c6d2917546..6fa7d05e37 100644 --- a/arc/mapping/engine_test.py +++ b/arc/mapping/engine_test.py @@ -57,7 +57,6 @@ def setUpClass(cls): (1.0669051052331406, -0.17519582095514982, 0.05416492980439295), (-0.6853171627400634, -0.8375353626879753, -0.028085652887100996), (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608))} - p_1_xyz = {'symbols': ('C', 'C', 'O', 'O', 'C', 'C', 'N', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 16, 12, 12, 14, 1, 1, 1, 1, 1, 1, 1, 1), 'coords': ((-2.4126223050522957, 1.1184165035580413, -0.4319265772396006), @@ -997,25 +996,24 @@ def test_flip_map(self): def test_make_bond_changes(self): """Test the make_bond_changes() function""" - spc1 = ARCSpecies(label="Test_bc", smiles="[CH2][CH2]") - spc2 = ARCSpecies(label="Test_bc", smiles="C=C") - label_species_atoms([spc1]) - label_species_atoms([spc2]) - rxn = ARCReaction(r_species = [spc1], p_species=[spc2]) - r_label_dict = {'*1': 0, '*2': 1} - l = [spc1] - make_bond_changes(rxn, l, r_label_dict) - self.assertTrue(spc2.mol.is_isomorphic(l[0].mol)) + # spc1 = ARCSpecies(label="Test_bc", smiles="[CH2][CH2]") + # spc2 = ARCSpecies(label="Test_bc", smiles="C=C") + # label_species_atoms([spc1]) + # label_species_atoms([spc2]) + # rxn = ARCReaction(r_species = [spc1], p_species=[spc2]) + # r_label_dict = {'*1': 0, '*2': 1} + # l = [spc1] + # make_bond_changes(rxn, l, r_label_dict) + # self.assertTrue(spc2.mol.is_isomorphic(l[0].mol)) + rxn = ARCReaction(r_species=[ARCSpecies(label="N4", smiles="NNNN")], p_species=[ARCSpecies(label="NH3", smiles="N"), ARCSpecies(label="NH2NHN_p", smiles="[N-]=[NH+]N")]) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn, backend="arc") - r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn, - rmg_reaction=rmg_reactions[0]) - assign_labels_to_products(rxn, p_label_dict) + rxn.family = '1,2_NH3_elimination' reactants, products = copy_species_list_for_mapping(rxn.r_species), copy_species_list_for_mapping(rxn.p_species) label_species_atoms(reactants), label_species_atoms(products) - r_bdes, _ = find_all_breaking_bonds(rxn, r_label_dict, True), find_all_breaking_bonds(rxn, p_label_dict, False) + print(f'\n\nfamily: {rxn.family}') + r_bdes, _ = find_all_breaking_bonds(rxn, True), find_all_breaking_bonds(rxn, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) self.assertFalse(r_cuts[1].mol.is_isomorphic(rxn.p_species[1].mol)) make_bond_changes(rxn=rxn, @@ -1052,14 +1050,12 @@ def test_update_xyz(self): def test_r_cut_p_cut_isomorphic(self): """Test the r_cut_p_cut_isomorphic() function""" rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) + reactants = copy_species_list_for_mapping(rxn_1_test.r_species) + products = copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, self.r_label_dict_rxn_1, True), find_all_breaking_bonds(rxn_1_test, self.p_label_dict_rxn_1, False) - + r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, True), find_all_breaking_bonds(rxn_1_test, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) - self.assertTrue(r_cut_p_cut_isomorphic(self.spc1,self.spc2)) for r_cut in r_cuts: for p_cut in p_cuts: @@ -1072,18 +1068,15 @@ def test_r_cut_p_cut_isomorphic(self): def test_pairing_reactants_and_products_for_mapping(self): """Test the pairing_reactants_and_products_for_mapping() function""" - smiles = ["[F]", "C[CH]C", "[CH3]"] + smiles = ['CC(=O)O[CH]CN', '[H]', '[CH3]'] rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) + reactants = copy_species_list_for_mapping(rxn_1_test.r_species) + products = copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, self.r_label_dict_rxn_1, True), find_all_breaking_bonds(rxn_1_test, self.p_label_dict_rxn_1, False) - + r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, True), find_all_breaking_bonds(rxn_1_test, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) - pairs_of_reactant_and_products = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts) - for pair in pairs_of_reactant_and_products: self.assertTrue(pair[0].mol.copy(deep=True).is_isomorphic(pair[1].mol.copy(deep=True))) self.assertIn(str(pair[0].mol.copy(deep=True).to_smiles()), smiles) @@ -1092,14 +1085,10 @@ def test_pairing_reactants_and_products_for_mapping(self): def test_map_pairs(self): """Test the map_pairs() function""" rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_1_test, backend="ARC") - r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_1_test, - rmg_reaction=rmg_reactions[0]) - reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) + reactants = copy_species_list_for_mapping(rxn_1_test.r_species) + products = copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, r_label_dict, True), find_all_breaking_bonds(rxn_1_test, p_label_dict, False) - + r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, True), find_all_breaking_bonds(rxn_1_test, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) pairs_of_reactant_and_products = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts) @@ -1116,28 +1105,26 @@ def test_map_pairs(self): self.assertEqual(map_[:3], [0, 1, 2]) self.assertIn(tuple(map_[3:6]), list(itertools.permutations([3, 4, 5]))) self.assertEqual(map_[6], 6) - self.assertIn(tuple(map_[7:]), list(itertools.permutations([7, 8, 9]))) + self.assertIn(tuple(map_[7:]), list(itertools.permutations([7, 8, 9, 10, 11, 12, 13, 14]))) self.assertEqual(len(np.unique(map_)), len(map_)) def test_glue_maps(self): + """Test the glue_maps() function""" rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_1_test, backend="ARC") - r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_1_test, - rmg_reaction=rmg_reactions[0]) - assign_labels_to_products(rxn_1_test, p_label_dict) - reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) + reactants = copy_species_list_for_mapping(rxn_1_test.r_species) + products = copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, r_label_dict, True), find_all_breaking_bonds(rxn_1_test, p_label_dict, False) - + r_bdes, p_bdes = find_all_breaking_bonds(rxn_1_test, True), find_all_breaking_bonds(rxn_1_test, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) pairs_of_reactant_and_products = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts) maps = map_pairs(pairs_of_reactant_and_products) - atom_map = glue_maps(maps,pairs_of_reactant_and_products) + atom_map = glue_maps(maps, pairs_of_reactant_and_products) self.assertEqual(len(atom_map), self.r_1.mol.get_num_atoms() + self.r_2.mol.get_num_atoms()) - atoms_r = [atom for atom in self.r_1.mol.copy(deep=True).atoms] + [atom for atom in self.r_2.mol.copy(deep=True).atoms] - atoms_p = [atom for atom in self.p_1.mol.copy(deep=True).atoms] + [atom for atom in self.p_2.mol.copy(deep=True).atoms] + atoms_r = [atom for atom in self.r_1.mol.copy(deep=True).atoms] + [atom for atom in + self.r_2.mol.copy(deep=True).atoms] + atoms_p = [atom for atom in self.p_1.mol.copy(deep=True).atoms] + [atom for atom in + self.p_2.mol.copy(deep=True).atoms] for index, value in enumerate(atom_map): self.assertEqual(atoms_r[index].symbol, atoms_p[value].symbol) @@ -1145,33 +1132,24 @@ def test_determine_bdes_on_spc_based_on_atom_labels(self): """tests the determine_bdes_indices_based_on_atom_labels() function""" spc = ARCSpecies(label="ethane", smiles="CC") label_species_atoms([spc]) - self.assertTrue(determine_bdes_on_spc_based_on_atom_labels(spc, (1, 2))) - self.assertFalse(determine_bdes_on_spc_based_on_atom_labels(spc, (2, 3))) + self.assertTrue(determine_bdes_on_spc_based_on_atom_labels(spc, (0, 1))) + self.assertFalse(determine_bdes_on_spc_based_on_atom_labels(spc, (1, 2))) def test_cut_species_based_on_atom_indices(self): """test the cut_species_based_on_atom_indices() function""" reactants = copy_species_list_for_mapping(self.arc_reaction_2.r_species) products = copy_species_list_for_mapping(self.arc_reaction_2.p_species) label_species_atoms(reactants), label_species_atoms(products) - product_dict = get_reaction_family_products(self.arc_reaction_2) - print(product_dict) - product_dict = product_dict[0] - r_label_dict_rxn_1 = product_dict['r_label_map'] - p_label_dict_rxn_1 = product_dict['p_label_map'] - print(r_label_dict_rxn_1) - print(p_label_dict_rxn_1) - - r_bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, label_dict=r_label_dict_rxn_1, r_direction=True) - p_bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, label_dict=p_label_dict_rxn_1, r_direction=False) - + # product_dicts = get_reaction_family_products(self.arc_reaction_2) + # product_dict = product_dicts[0] + r_bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, r_direction=True) + print(f'product atoms of rxn2: 0: {self.arc_reaction_2.p_species[0].mol.atoms}, 1: {self.arc_reaction_2.p_species[1].mol.atoms}') + p_bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, r_direction=False) + print(f'r_bdes: {r_bdes}, p_bdes: {p_bdes}') + print(f'p atoms: 0: {products[0].mol.atoms}, 1: {products[1].mol.atoms}') r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) - - print(f'r_bdes: {r_bdes}') - print(f'p_bdes: {p_bdes}') - print(f'r_cuts: {[s.mol.to_smiles() for s in r_cuts]}') - print(f'p_cuts: {[s.mol.to_smiles() for s in p_cuts]}') - + print(f'r_cuts: {[r.copy().mol.to_smiles() for r in r_cuts]}, p_cuts: {[r.copy().mol.to_smiles() for r in p_cuts]}') self.assertIn("[CH2]CC", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) self.assertIn("[H]", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) self.assertIn("[NH2]", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) @@ -1179,10 +1157,10 @@ def test_cut_species_based_on_atom_indices(self): self.assertIn("[H]", [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) self.assertIn("[NH2]", [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) - spc = ARCSpecies(label="test", smiles="CNC", bdes = [(1, 2), (2, 3)]) + spc = ARCSpecies(label="test", smiles="CNC") for i, a in enumerate(spc.mol.atoms): a.label=str(i) - cuts = cut_species_based_on_atom_indices([spc], [(1, 2), (2, 3)]) + cuts = cut_species_based_on_atom_indices([spc], [(0, 1), (1, 2)]) self.assertEqual(len(cuts), 3) for cut in cuts: self.assertTrue(any([cut.mol.copy(deep=True).is_isomorphic(ARCSpecies(label="1", smiles="[CH3]").mol), @@ -1191,34 +1169,94 @@ def test_cut_species_based_on_atom_indices(self): h2 = ARCSpecies(label="H2", smiles="[H][H]") label_species_atoms([h2]) - cuts = cut_species_based_on_atom_indices([h2], [(1, 2)]) + cuts = cut_species_based_on_atom_indices([h2], [(0, 1)]) self.assertEqual(len(cuts), 2) for cut in cuts: self.assertEqual(cut.get_xyz()["symbols"], ('H',)) - spcs = [ARCSpecies(label="r", smiles = 'O=C(O)CCF')] + spcs = [ARCSpecies(label="r", smiles = 'O=C(O)CCF', + xyz={'symbols': ('O', 'C', 'O', 'C', 'C', 'F', 'H', 'H', 'H', 'H', 'H'), + 'isotopes': (16, 12, 16, 12, 12, 19, 1, 1, 1, 1, 1), + 'coords': ((1.1567401394807186, 0.3292352797648203, -1.4438840497688), + (1.079794258193769, 0.18837439525248412, -0.23580436209026492), + (2.184999596455945, -0.02862572694090521, 0.5028447307908376), + (-0.15828046358747527, 0.22516852199733842, 0.6026621014486843), + (-1.3729926295596129, 0.4703126054291076, -0.26280242788844443), + (-2.4751613391679754, 0.4966265394383556, 0.5356539921809551), + (2.9199184814309573, -0.03170785916941691, -0.14635137063038064), + (-0.26755118310691867, -0.7279235396445063, 1.1300570091365727), + (-0.0686732659284988, 1.0230262342082468, 1.3468825326130487), + (-1.313896668558704, 1.4294967806185437, -0.7865558706137695), + (-1.5129630024826421, -0.32309516060704263, -1.0035834767256502))})] label_species_atoms(spcs) - cuts = cut_species_based_on_atom_indices(spcs, [(6, 5), (4, 2), (3, 7)]) + cuts = cut_species_based_on_atom_indices(spcs, [(1, 2), (3, 4), (4, 10)]) self.assertEqual(len(cuts), 4) - def test_translate_indices_based_on_ref_species(self): - """tests the translate_indices_based_on_ref_species() function""" + spcs = [ARCSpecies(label="ring", smiles = 'CC1CCOCC1[N+](=O)([O-])', + xyz={'symbols': ('C', 'C', 'C', 'C', 'O', 'C', 'C', 'N', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H'), + 'isotopes': (12, 12, 12, 12, 16, 12, 12, 14, 16, 16, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), + 'coords': ((-0.19963202410071415, 2.11431236883385, 0.35512152305302197), + (-0.21499972296989353, 0.5855566636041595, 0.4252779409474618), + (-1.5279020878978917, 0.042997223439599976, -0.16475709380965647), + (-1.5142756264743436, -1.4801435389166957, -0.23348399770872003), + (-0.39427404620551165, -1.9345051360499914, -0.9947303725939133), + (0.8504594323376896, -1.5617817457794443, -0.3915176374284775), + (0.9727667641076266, -0.04585137142160857, -0.3203888579911006), + (2.266478114387414, 0.2730161804968195, 0.3742376225644008), + (2.395228629259047, -0.08871684931314952, 1.5513228364332547), + (3.1195221496216616, 0.8781087874071717, -0.2865464895135631), + (-1.0748919586201595, 2.532876649227796, 0.8634243187970725), + (-0.2125798569006862, 2.465012686844741, -0.6824090093235058), + (0.6909911289491314, 2.5252788705772518, 0.8408874242217174), + (-0.18391181140032295, 0.2986578719375753, 1.4853677298816432), + (-1.6551845431159757, 0.44049461022166125, -1.1800168508998137), + (-2.3833437995959468, 0.3815988059611917, 0.4305386229780667), + (-1.4830042838432955, -1.925925263188267, 0.7675312553931687), + (-2.4218853597129604, -1.8388555520124514, -0.7290646661080117), + (0.9188046169404014, -2.040202832179712, 0.5936529337480443), + (1.6528379049665831, -1.9828525321277448, -1.0079109263010095), + (1.0337508478988304, 0.36617104486181684, -1.3359287392782864))})] + label_species_atoms(spcs) + cuts = cut_species_based_on_atom_indices(spcs, [(1, 2), (4, 5)]) + self.assertEqual(len(cuts), 2) + + def test_translate_bdes_based_on_ref_species(self): + """tests the translate_bdes_based_on_ref_species() function""" spc_1 = ARCSpecies(label='pentanol', smiles='CCCO') spc_2 = spc_1.copy() - new_bdes = translate_indices_based_on_ref_species(species=[spc_1], ref_species=[spc_2], bdes=[(0, 1)]) + new_bdes = translate_bdes_based_on_ref_species(species=[spc_1], ref_species=[spc_2], bdes=[(0, 1)]) self.assertEqual(new_bdes, [(0, 1)]) spc_1 = ARCSpecies(label='C3Oj', smiles='CCC[O]') spc_2 = ARCSpecies(label='H', smiles='[H]') spc_3 = ARCSpecies(label='OH', smiles='[OH]') - new_bdes = translate_indices_based_on_ref_species(species=[spc_3, spc_2, spc_1], - ref_species=[spc_3, spc_2, spc_1], - bdes=[(0, 2), (2, 6)]) + new_bdes = translate_bdes_based_on_ref_species(species=[spc_3, spc_2, spc_1], + ref_species=[spc_3, spc_2, spc_1], + bdes=[(0, 2), (2, 6)]) self.assertEqual(new_bdes, [(0, 2), (2, 6)]) - new_bdes = translate_indices_based_on_ref_species(species=[spc_1, spc_2, spc_3], - ref_species=[spc_3, spc_2, spc_1], - bdes=[(0, 2), (2, 6)]) - self.assertEqual(new_bdes, [(12, 11), (11, 3)]) + new_bdes = translate_bdes_based_on_ref_species(species=[spc_1, spc_2, spc_3], + ref_species=[spc_3, spc_2, spc_1], + bdes=[(0, 2), (2, 6)]) + self.assertEqual(new_bdes, [(11, 12), (3, 11)]) + + def test_translate_indices_based_on_ref_species(self): + """tests the translate_index_based_on_ref_species() function""" + spc_1 = ARCSpecies(label='pentanol', smiles='CCCO') + spc_2 = spc_1.copy() + new_indices = translate_indices_based_on_ref_species(species=[spc_1], ref_species=[spc_2], indices=[0]) + self.assertEqual(new_indices, [0]) + + spc_1 = ARCSpecies(label='C3Oj', smiles='CCC[O]') + spc_2 = ARCSpecies(label='H', smiles='[H]') + spc_3 = ARCSpecies(label='OH', smiles='[OH]') + new_indices = translate_indices_based_on_ref_species(species=[spc_3, spc_2, spc_1], + ref_species=[spc_3, spc_2, spc_1], + indices=[0, 1, 2, 10]) + self.assertEqual(new_indices, [0, 1, 2, 10]) + new_indices = translate_indices_based_on_ref_species(species=[spc_1, spc_2, spc_3], + ref_species=[spc_3, spc_2, spc_1], + indices=[0, 1, 2, 10]) + self.assertEqual(new_indices, [12, 13, 11, 7]) def test_copy_species_list_for_mapping(self): """tests the copy_species_list_for_mapping() function""" @@ -1235,25 +1273,23 @@ def test_copy_species_list_for_mapping(self): def test_find_all_breaking_bonds(self): """tests the find_all_breaking_bonds() function""" rxn = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - products_dicts = get_reaction_family_products(rxn) - bdes = find_all_breaking_bonds(rxn=rxn, label_dict=products_dicts[0]['r_label_map'], r_direction=True) + bdes = find_all_breaking_bonds(rxn=rxn, r_direction=True) self.assertEqual(bdes, [(4, 10)]) - bdes = find_all_breaking_bonds(rxn=rxn, label_dict=products_dicts[0]['r_label_map'], r_direction=False) - self.assertEqual(bdes, [(10, 16)]) + bdes = find_all_breaking_bonds(rxn=rxn, r_direction=False) + self.assertEqual(bdes, [(15, 19)]) - products_dicts = get_reaction_family_products(self.arc_reaction_1) - bdes = find_all_breaking_bonds(rxn=self.arc_reaction_1, label_dict=products_dicts[0]['r_label_map'], r_direction=True) + bdes = find_all_breaking_bonds(rxn=self.arc_reaction_1, r_direction=True) self.assertEqual(bdes, [(0, 1)]) - products_dicts = get_reaction_family_products(self.arc_reaction_2) - bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, label_dict=products_dicts[0]['r_label_map'], r_direction=True) + bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, r_direction=True) self.assertEqual(bdes, [(0, 3)]) - bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, label_dict=products_dicts[0]['r_label_map'], r_direction=False) - self.assertEqual(bdes, [(3, 11)]) + bdes = find_all_breaking_bonds(rxn=self.arc_reaction_2, r_direction=False) + print(f'p atoms: 0: {self.arc_reaction_2.p_species[0].mol.atoms}, 1: {self.arc_reaction_2.p_species[1].mol.atoms}') # todo: fix product BDE + self.assertEqual(bdes, [(10, 13)]) + raise - products_dicts = get_reaction_family_products(self.arc_reaction_3) - bdes = find_all_breaking_bonds(rxn=self.arc_reaction_3, label_dict=products_dicts[0]['r_label_map'], r_direction=False) - self.assertEqual(bdes, [(0, 4)]) + bdes = find_all_breaking_bonds(rxn=self.arc_reaction_3, r_direction=False) + self.assertEqual(bdes, [(0, 1)]) if __name__ == '__main__': diff --git a/arc/processor.py b/arc/processor.py index 39ca0b74d3..afcc7851d8 100644 --- a/arc/processor.py +++ b/arc/processor.py @@ -6,10 +6,9 @@ import shutil from typing import Optional, Type -from rmgpy.data.rmg import RMGDatabase +from rmgpy.data.rmg import RMGDatabase # todo import arc.plotter as plotter -import arc.rmgdb as rmgdb from arc.common import get_logger from arc.level import Level from arc.statmech.factory import statmech_factory @@ -265,7 +264,7 @@ def compare_rates(rxns_for_kinetics_lib: list, reactions_to_compare = list() # reactions for which a rate was both calculated and estimated. for reaction in rxns_for_kinetics_lib: try: - reaction.rmg_reactions = rmgdb.determine_rmg_kinetics(rmgdb=rmg_database, + reaction.rmg_reactions = rmgdb.determine_rmg_kinetics(rmgdb=rmg_database, # todo reaction=reaction.rmg_reaction, dh_rxn298=reaction.dh_rxn298) except Exception as e: @@ -317,7 +316,7 @@ def load_rmg_database(rmg_database: Optional[Type[RMGDatabase]], for species in species_dict.values()]): load_thermo_libs = True if rmg_database is not None and (load_kinetic_libs or load_thermo_libs): - rmgdb.load_rmg_database(rmgdb=rmg_database, + rmgdb.load_rmg_database(rmgdb=rmg_database, # todo load_thermo_libs=load_thermo_libs, load_kinetic_libs=load_kinetic_libs) diff --git a/arc/reaction/reaction.py b/arc/reaction/reaction.py index ab53320853..4042cf5470 100644 --- a/arc/reaction/reaction.py +++ b/arc/reaction/reaction.py @@ -581,7 +581,7 @@ def determine_family(self, consider_rmg_families: bool = True, consider_arc_families: bool = True, discover_own_reverse_rxns_in_reverse: bool = False, - ): + ) -> Tuple[Optional[str], Optional[bool]]: """ Determine the RMG reaction family. Populates the .family, and .family_own_reverse attributes. @@ -591,6 +591,9 @@ def determine_family(self, consider_rmg_families (bool, optional): Whether to consider RMG's families in addition to ARC's. consider_arc_families (bool, optional): Whether to consider ARC's families in addition to RMG's. discover_own_reverse_rxns_in_reverse (bool, optional): Whether to discover own reverse reactions in reverse. + + Returns: + Tuple[Optional[str], Optional[bool]]: The RMG reaction family and the family's own_reverse property. """ if self.rmg_reaction is not None: product_dicts = get_reaction_family_products(rxn=self, diff --git a/arc/species/converter.py b/arc/species/converter.py index a813381f38..1aaaa79b06 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -1899,6 +1899,25 @@ def check_isomorphism(mol1, mol2, filter_structures=True, convert_to_single_bond return False +def check_molecule_list_order(mols_1: List[Molecule], mols_2: List[Molecule]): + """ + Check if the order of molecules in two lists is the same. + + Args: + mols_1 (List[Molecule]): A list of RMG Molecule objects. + mols_2 (List[Molecule]): A list of RMG Molecule objects. + + Returns: + bool: Whether the order of molecules in the two lists is the same. + """ + if len(mols_1) != len(mols_2): + return False + for mol1, mol2 in zip(mols_1, mols_2): + if not check_isomorphism(mol1, mol2): + return False + return True + + def get_center_of_mass(xyz): """ Get the center of mass of xyz coordinates. diff --git a/arc/species/species_test.py b/arc/species/species_test.py index 2813e840d4..69eaeb119f 100644 --- a/arc/species/species_test.py +++ b/arc/species/species_test.py @@ -2599,7 +2599,7 @@ def test_correct_representation_of_o2_and_s2(self): o2 = ARCSpecies(label='O2', smiles='[O][O]', xyz="""O 0.0000000 0.0000000 0.6029240 O 0.0000000 0.0000000 -0.6029240""") for mol in o2.mol_list: - print(f'mol lost') + print(f'mol list') print(mol.to_smiles()) self.assertEqual(o2.multiplicity, 3) self.assertEqual(o2.mol.to_smiles(), '[O][O]') @@ -2642,6 +2642,7 @@ def test_split_mol(self): self.assertEqual(m.to_smiles(), 'O') self.assertEqual(fragments, [[0, 3, 4], [1, 5, 6], [2, 7, 8]]) + @classmethod def tearDownClass(cls): """ diff --git a/arc/testing/test_JobAdapter_ServerTimeLimit/calcs/Species/spc1/spc1/submit.sh b/arc/testing/test_JobAdapter_ServerTimeLimit/calcs/Species/spc1/spc1/submit.sh new file mode 100644 index 0000000000..658bc41607 --- /dev/null +++ b/arc/testing/test_JobAdapter_ServerTimeLimit/calcs/Species/spc1/spc1/submit.sh @@ -0,0 +1,37 @@ +#!/bin/bash -l +#SBATCH -p normal +#SBATCH -J server3 +#SBATCH -N 1 +#SBATCH -n 8 +#SBATCH --time=120:00:00 +#SBATCH --mem-per-cpu=15770000000 +#SBATCH -o out.txt +#SBATCH -e err.txt + +export g16root=/home/gridsan/groups/GRPAPI/Software +export PATH=$g16root/g16/:$g16root/gv:$PATH +which g16 + +echo "============================================================" +echo "Job ID : $SLURM_JOB_ID" +echo "Job Name : $SLURM_JOB_NAME" +echo "Starting on : $(date)" +echo "Running on node : $SLURMD_NODENAME" +echo "Current directory : $(pwd)" +echo "============================================================" + +touch initial_time + +GAUSS_SCRDIR=/state/partition1/user//$SLURM_JOB_NAME-$SLURM_JOB_ID +export $GAUSS_SCRDIR +. $g16root/g16/bsd/g16.profile + +mkdir -p $GAUSS_SCRDIR + +g16 < input.gjf > input.log + +rm -rf $GAUSS_SCRDIR + +touch final_time + + \ No newline at end of file