diff --git a/arc/job/adapters/ts/gcn_ts.py b/arc/job/adapters/ts/gcn_ts.py index 1e2015f6c1..0e549cfe4e 100644 --- a/arc/job/adapters/ts/gcn_ts.py +++ b/arc/job/adapters/ts/gcn_ts.py @@ -274,46 +274,53 @@ def execute_gcn(self, exe_type: str = 'incore'): charge=rxn.charge, multiplicity=rxn.multiplicity, ) - write_sdf_files(rxn=rxn, - reactant_path=self.reactant_path, - product_path=self.product_path, - ) - if exe_type == 'queue': - input_dict = {'reactant_path': self.reactant_path, - 'product_path': self.product_path, - 'local_path': self.local_path, - 'yml_out_path': self.yml_out_path, - 'repetitions': self.repetitions, - } - save_yaml_file(path=self.yml_in_path, content=input_dict) - self.legacy_queue_execution() - elif exe_type == 'incore': - for _ in range(self.repetitions): - run_subprocess_locally(direction='F', - reactant_path=self.reactant_path, - product_path=self.product_path, - ts_path=self.ts_fwd_path, - local_path=self.local_path, - ts_species=rxn.ts_species, - ) - run_subprocess_locally(direction='R', - reactant_path=self.product_path, - product_path=self.reactant_path, - ts_path=self.ts_rev_path, - local_path=self.local_path, - ts_species=rxn.ts_species, - ) - if len(self.reactions) < 5: - successes = len([tsg for tsg in rxn.ts_species.ts_guesses if tsg.success and 'gcn' in tsg.method]) - if successes: - logger.info(f'GCN successfully found {successes} TS guesses for {rxn.label}.') - else: - logger.info(f'GCN did not find any successful TS guesses for {rxn.label}.') + for i in range(len(rxn.atom_maps)): + try: + write_sdf_files(rxn=rxn, + reactant_path=self.reactant_path, + product_path=self.product_path, + am_index=i, + ) + except IndexError: + logger.warning(f'GCN adapter could not write SDF files for {rxn.label} with atom map {i}.') + continue + if exe_type == 'queue': + input_dict = {'reactant_path': self.reactant_path, + 'product_path': self.product_path, + 'local_path': self.local_path, + 'yml_out_path': self.yml_out_path, + 'repetitions': self.repetitions, + } + save_yaml_file(path=self.yml_in_path, content=input_dict) + self.legacy_queue_execution() + elif exe_type == 'incore': + for _ in range(self.repetitions): + run_subprocess_locally(direction='F', + reactant_path=self.reactant_path, + product_path=self.product_path, + ts_path=self.ts_fwd_path, + local_path=self.local_path, + ts_species=rxn.ts_species, + ) + run_subprocess_locally(direction='R', + reactant_path=self.product_path, + product_path=self.reactant_path, + ts_path=self.ts_rev_path, + local_path=self.local_path, + ts_species=rxn.ts_species, + ) + if len(self.reactions) < 5: + successes = len([tsg for tsg in rxn.ts_species.ts_guesses if tsg.success and 'gcn' in tsg.method]) + if successes: + logger.info(f'GCN successfully found {successes} TS guesses for {rxn.label}.') + else: + logger.info(f'GCN did not find any successful TS guesses for {rxn.label}.') def write_sdf_files(rxn: 'ARCReaction', reactant_path: str, product_path: str, + am_index: int = 0, ): """ Write reactant and product SDF files using RDKit. @@ -322,9 +329,10 @@ def write_sdf_files(rxn: 'ARCReaction', rxn (ARCReaction): The relevant reaction. reactant_path (str): The path to the reactant SDF file. product_path (str): The path to the product SDF file. + am_index (int, optional): The atom map index. Default: 0. """ reactant_rdkit_mol = rdkit_conf_from_mol(rxn.r_species[0].mol, rxn.r_species[0].get_xyz())[1] - mapped_product = rxn.get_single_mapped_product_xyz() + mapped_product = rxn.get_single_mapped_product_xyz(am_index) product_rdkit_mol = rdkit_conf_from_mol(mapped_product.mol, mapped_product.get_xyz())[1] w = Chem.SDWriter(reactant_path) w.write(reactant_rdkit_mol) diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py index 19e1aca401..cf77ff77b7 100644 --- a/arc/job/adapters/ts/heuristics_test.py +++ b/arc/job/adapters/ts/heuristics_test.py @@ -302,10 +302,10 @@ def test_heuristics_for_h_abstraction(self): products=[Species(smiles='[CH3]'), Species(smiles='[H][H]')])) rxn4.determine_family(rmg_database=self.rmgdb) self.assertEqual(rxn4.family.label, 'H_Abstraction') - self.assertEqual(rxn4.atom_map[0], 0) + self.assertEqual(rxn4.atom_maps[0][0], 0) for index in [1, 2, 3, 4]: - self.assertIn(rxn4.atom_map[index], [1, 2, 3, 4, 5]) - self.assertIn(rxn4.atom_map[5], [4, 5]) + self.assertIn(rxn4.atom_maps[0][index], [1, 2, 3, 4, 5]) + self.assertIn(rxn4.atom_maps[0][5], [4, 5]) heuristics_4 = HeuristicsAdapter(job_type='tsg', reactions=[rxn4], testing=True, @@ -926,13 +926,13 @@ def test_keeping_atom_order_in_ts(self): p_species=[ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz), ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz)]) rxn_1.determine_family(rmg_database=self.rmgdb) - self.assertIn(rxn_1.atom_map[0], [0, 1]) - self.assertIn(rxn_1.atom_map[1], [0, 1]) + self.assertIn(rxn_1.atom_maps[0][0], [0, 1]) + self.assertIn(rxn_1.atom_maps[0][1], [0, 1]) for index in [2, 3, 4, 5, 6, 7]: - self.assertIn(rxn_1.atom_map[index], [2, 3, 4, 5, 6, 16]) - self.assertEqual(rxn_1.atom_map[8:12], [7, 8, 9, 10]) - self.assertIn(tuple(rxn_1.atom_map[12:15]), itertools.permutations([13, 11, 12])) - self.assertIn(rxn_1.atom_map[15:], [[14, 15], [15, 14]]) + self.assertIn(rxn_1.atom_maps[0][index], [2, 3, 4, 5, 6, 16]) + self.assertEqual(rxn_1.atom_maps[0][8:12], [7, 8, 9, 10]) + self.assertIn(tuple(rxn_1.atom_maps[0][12:15]), itertools.permutations([13, 11, 12])) + self.assertIn(rxn_1.atom_maps[0][15:], [[14, 15], [15, 14]]) heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn_1], testing=True, @@ -952,12 +952,12 @@ def test_keeping_atom_order_in_ts(self): ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz)]) rxn_2.determine_family(rmg_database=self.rmgdb) self.assertEqual(rxn_2.family.label, 'H_Abstraction') - self.assertEqual(rxn_2.atom_map[:2], [11, 10]) - self.assertIn(tuple(rxn_2.atom_map[2:5]), itertools.permutations([9, 16, 15])) - self.assertIn(tuple(rxn_2.atom_map[5:8]), itertools.permutations([12, 13, 14])) - self.assertEqual(rxn_2.atom_map[8:12], [0, 1, 2, 3]) - self.assertIn(tuple(rxn_2.atom_map[12:15]), itertools.permutations([4, 5, 6])) - self.assertIn(tuple(rxn_2.atom_map[15:]), itertools.permutations([7, 8])) + self.assertEqual(rxn_2.atom_maps[0][:2], [11, 10]) + self.assertIn(tuple(rxn_2.atom_maps[0][2:5]), itertools.permutations([9, 16, 15])) + self.assertIn(tuple(rxn_2.atom_maps[0][5:8]), itertools.permutations([12, 13, 14])) + self.assertEqual(rxn_2.atom_maps[0][8:12], [0, 1, 2, 3]) + self.assertIn(tuple(rxn_2.atom_maps[0][12:15]), itertools.permutations([4, 5, 6])) + self.assertIn(tuple(rxn_2.atom_maps[0][15:]), itertools.permutations([7, 8])) heuristics_2 = HeuristicsAdapter(job_type='tsg', reactions=[rxn_2], testing=True, @@ -976,12 +976,12 @@ def test_keeping_atom_order_in_ts(self): p_species=[ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz), ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz)]) rxn_3.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn_3.atom_map[:4], [7, 8, 9, 10]) - self.assertIn(tuple(rxn_3.atom_map[4:7]), itertools.permutations([11, 12, 13])) - self.assertIn(tuple(rxn_3.atom_map[7:9]), itertools.permutations([14, 15])) - self.assertEqual(rxn_3.atom_map[9:11], [1, 0]) - self.assertIn(tuple(rxn_3.atom_map[11:14]), itertools.permutations([16, 5, 6])) - self.assertIn(tuple(rxn_3.atom_map[14:]), itertools.permutations([3, 4, 2])) + self.assertEqual(rxn_3.atom_maps[0][:4], [7, 8, 9, 10]) + self.assertIn(tuple(rxn_3.atom_maps[0][4:7]), itertools.permutations([11, 12, 13])) + self.assertIn(tuple(rxn_3.atom_maps[0][7:9]), itertools.permutations([14, 15])) + self.assertEqual(rxn_3.atom_maps[0][9:11], [1, 0]) + self.assertIn(tuple(rxn_3.atom_maps[0][11:14]), itertools.permutations([16, 5, 6])) + self.assertIn(tuple(rxn_3.atom_maps[0][14:]), itertools.permutations([3, 4, 2])) heuristics_3 = HeuristicsAdapter(job_type='tsg', reactions=[rxn_3], @@ -1001,12 +1001,12 @@ def test_keeping_atom_order_in_ts(self): p_species=[ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz), ARCSpecies(label='C2H5', smiles='C[CH2]', xyz=self.c2h5_xyz)]) rxn_4.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn_4.atom_map[:4], [0, 1, 2, 3]) - self.assertIn(tuple(rxn_4.atom_map[4:7]), itertools.permutations([4, 5, 6])) - self.assertIn(tuple(rxn_4.atom_map[7:9]), itertools.permutations([7, 8])) - self.assertEqual(rxn_4.atom_map[9:11], [11, 10]) - self.assertIn(tuple(rxn_4.atom_map[11:14]), itertools.permutations([9, 15, 16])) - self.assertIn(tuple(rxn_4.atom_map[14:]), itertools.permutations([12, 13, 14 ])) + self.assertEqual(rxn_4.atom_maps[0][:4], [0, 1, 2, 3]) + self.assertIn(tuple(rxn_4.atom_maps[0][4:7]), itertools.permutations([4, 5, 6])) + self.assertIn(tuple(rxn_4.atom_maps[0][7:9]), itertools.permutations([7, 8])) + self.assertEqual(rxn_4.atom_maps[0][9:11], [11, 10]) + self.assertIn(tuple(rxn_4.atom_maps[0][11:14]), itertools.permutations([9, 15, 16])) + self.assertIn(tuple(rxn_4.atom_maps[0][14:]), itertools.permutations([12, 13, 14 ])) heuristics_4 = HeuristicsAdapter(job_type='tsg', reactions=[rxn_4], testing=True, diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index e378eedd69..f7e0caa91e 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -40,64 +40,63 @@ def map_reaction(rxn: 'ARCReaction', backend: str = 'ARC', db: Optional['RMGDatabase'] = None, - flip = False - ) -> Optional[List[int]]: + flip: bool = False + ) -> Optional[List[List[int]]]: """ - Map a reaction. + Map a reaction. A wrapper for the ``map_rxn`` function. 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. + flip (bool, optional): Try mapping with a flipped reaction. Returns: - Optional[List[int]]: + Optional[List[List[int]]]: Entry indices are running atom indices of the reactants, corresponding entry values are running atom indices of the products. """ if flip: - logger.warning(f"The requested ARC reaction {rxn} could not be atom mapped using {backend}. Trying again with the flipped reaction.") + logger.warning(f"The requested ARC reaction {rxn} could not be atom mapped using {backend}. " + f"Trying again with the flipped reaction.") try: - _map = flip_map(map_rxn(rxn.flip_reaction(), backend=backend, db=db)) + atom_maps = flip_maps(map_rxn(rxn.flip_reaction(), backend=backend)) except ValueError: return None - return _map + return atom_maps else: if rxn.family is None: rmgdb.determine_family(reaction=rxn, db=db) if rxn.family is None: - logger.warning(f'Could not determine the reaction family for {rxn.label}. Mapping as a general or isomerization reaction.') - _map = map_general_rxn(rxn, backend=backend) - return _map if _map is not None else map_reaction(rxn, backend=backend, db=db, flip=True) + logger.warning(f'Could not determine the reaction family for {rxn.label}.\n' + f'Mapping as a general or isomerization reaction.') + atom_maps = map_general_rxn(rxn) + return atom_maps if atom_maps is not None else map_reaction(rxn, backend=backend, db=db, flip=True) try: - _map = map_rxn(rxn, backend=backend, db=db) - except ValueError as e: + atom_maps = map_rxn(rxn, backend=backend) + except ValueError: return map_reaction(rxn, backend=backend, db=db, flip=True) - return _map if _map is not None else map_reaction(rxn, backend=backend, db=db, flip=True) + return atom_maps if atom_maps is not None else map_reaction(rxn, backend=backend, db=db, flip=True) def map_general_rxn(rxn: 'ARCReaction', - backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, - ) -> Optional[List[int]]: + ) -> Optional[List[List[int]]]: """ - Map a general reaction (one that was not categorized into a reaction family by RMG). - The general method isn't great, a family-specific method should be implemented where possible. + Map a general reaction that could not be categorized into an RMG reaction family. + The general method isn't great, using the family-based method is better. 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]]: + Optional[List[List[int]]]: Entry indices are running atom indices of the reactants, corresponding entry values are running atom indices of the products. """ if rxn.is_isomerization(): - atom_map = map_isomerization_reaction(rxn=rxn) - if atom_map is not None: - return atom_map + atom_maps = map_isomerization_reaction(rxn=rxn) + if atom_maps is not None: + return atom_maps # If the reaction is not a known RMG template and is not isomerization, use fragments via the QCElemental backend. qcmol_1 = create_qc_mol(species=[spc.copy() for spc in rxn.r_species], @@ -112,21 +111,19 @@ def map_general_rxn(rxn: 'ARCReaction', return None data = qcmol_2.align(ref_mol=qcmol_1, verbose=0)[1] atom_map = data['mill'].atommap.tolist() - return atom_map + return [atom_map] def map_isomerization_reaction(rxn: 'ARCReaction', - ) -> Optional[List[int]]: + ) -> Optional[List[List[int]]]: """ Map isomerization reaction that has no corresponding RMG family. 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]]: + Optional[List[List[int]]]: Entry indices are running atom indices of the reactants, corresponding entry values are running atom indices of the products. """ @@ -203,18 +200,17 @@ def map_isomerization_reaction(rxn: 'ARCReaction', for bond in atom.bonds.values(): bond.order = 1 if success: - return map_two_species(r_copy, p_copy, map_type='list') + return [map_two_species(r_copy, p_copy, map_type='list')] # Fallback to the general mapping algorithm. - return map_two_species(rxn.r_species[0], rxn.p_species[0], map_type='list') + return [map_two_species(rxn.r_species[0], rxn.p_species[0], map_type='list')] def map_rxn(rxn: 'ARCReaction', backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, - ) -> Optional[List[int]]: + ) -> Optional[List[List[int]]]: """ - A wrapper function for mapping reaction, uses databases for mapping with the correct reaction family parameters. + Map a reaction. The main function of the atom mapping module. Strategy: 1) get_rmg_reactions_from_arc_reaction, get_atom_indices_of_labeled_atoms_in_an_rmg_reaction. 2) (For bimolecular reactions) Find the species in which the bond is broken. @@ -224,50 +220,47 @@ def map_rxn(rxn: 'ARCReaction', 6) Join maps together. Args: - rxn (ARCReaction): An ARCReaction object instance that belongs to the RMG H_Abstraction reaction family. + 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]]: + Optional[List[List[int]]]: Entry indices are running atom indices of the reactants, 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) - + if not rmg_reactions: return None - - r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn, - rmg_reaction=rmg_reactions[0]) - # step 2: - assign_labels_to_products(rxn, p_label_dict) - - #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) + atom_maps = list() + for rmg_reaction in rmg_reactions: + r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn, + rmg_reaction=rmg_reaction) - r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) - p_cuts = cut_species_based_on_atom_indices(products, p_bdes) + assign_labels_to_products(rxn, p_label_dict) + 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) - try: - make_bond_changes(rxn, r_cuts, r_label_dict) - except (ValueError, IndexError, ActionError, AtomTypeError) as e: - logger.warning(e) + r_bdes, p_bdes = find_all_bdes(rxn, r_label_dict, True), find_all_bdes(rxn, p_label_dict, False) + r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) + p_cuts = cut_species_based_on_atom_indices(products, p_bdes) - r_cuts, p_cuts = update_xyz(r_cuts), update_xyz(p_cuts) + try: + make_bond_changes(rxn, r_cuts, r_label_dict) + except (ValueError, IndexError, ActionError, AtomTypeError) as e: + logger.warning(e) + r_cuts, p_cuts = update_xyz(r_cuts), update_xyz(p_cuts) - #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]}") - return None - # step 5: - maps = map_pairs(pairs_of_reactant_and_products) + 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]}") + return None + + atom_map = map_pairs(pairs_of_reactant_and_products) + if atom_map is not None: + atom_map = glue_maps(atom_map, pairs_of_reactant_and_products) + atom_maps.append(atom_map) - #step 6: - return glue_maps(maps, pairs_of_reactant_and_products) + return atom_maps diff --git a/arc/mapping/engine.py b/arc/mapping/engine.py index 45d97ce724..b30c1d00b6 100644 --- a/arc/mapping/engine.py +++ b/arc/mapping/engine.py @@ -267,7 +267,7 @@ def map_two_species(spc_1: Union[ARCSpecies, Species, Molecule], return {0: 0} return [0] - # A shortcut for homonuclear diatomic species. + # A shortcut for homo-nuclear diatomic species. if spc_1.number_of_atoms == spc_2.number_of_atoms == 2 \ and len(set([atom.element.symbol for atom in spc_1.mol.atoms])) == 1: if map_type == 'dict': @@ -288,65 +288,105 @@ def map_two_species(spc_1: Union[ARCSpecies, Species, Molecule], if backend.lower() not in ['qcelemental', 'arc']: raise ValueError(f'The backend method could be either "QCElemental" or "ARC", got {backend}.') - atom_map = None + atom_map = None if backend.lower() == 'arc': + atom_map = map_two_species_via_arc(spc_1=spc_1, + spc_2=spc_2, + map_type=map_type, + consider_chirality=consider_chirality, + ) + if atom_map is None and allow_backend_shift: + backend = 'QCElemental' + if backend.lower() == 'qcelemental': + atom_map = map_two_species_via_qcel(spc_1=spc_1, spc_2=spc_2, map_type=map_type) + + if inc_vals is not None: + atom_map = [value + inc_vals for value in atom_map] + return atom_map + + +def map_two_species_via_arc(spc_1: Union[ARCSpecies, Species, Molecule], + spc_2: Union[ARCSpecies, Species, Molecule], + map_type: str = 'list', + consider_chirality: bool = True, + ) -> Optional[Union[List[int], Dict[int, int]]]: + """ + Map the atoms in spc_1 to the atoms in spc_2 using the ARC 3DAM method. + + Args: + spc_1 (Union[ARCSpecies, Species, Molecule]): Species 1. + spc_2 (Union[ARCSpecies, Species, Molecule]): Species 2. + map_type (str, optional): Whether to return a 'list' or a 'dict' map type. + consider_chirality (bool, optional): Whether to consider chirality when fingerprinting. + + Returns: + Optional[Union[List[int], Dict[int, int]]]: + The atom maps. By default, a list of atom maps is returned. + """ + fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) + fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) + candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) + if candidates is None or len(candidates) == 0: + consider_chirality = not consider_chirality fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) if candidates is None or len(candidates) == 0: - consider_chirality = not consider_chirality - fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) - fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) - candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) - if candidates is None or len(candidates) == 0: - logger.warning(f'Could not identify superimposable candidates {spc_1} and {spc_2}.') - return None - if not len(candidates): - if allow_backend_shift: - backend = 'QCElemental' - else: - return None - else: - rmsds, fixed_spcs = list(), list() - for candidate in candidates: - fixed_spc_1, fixed_spc_2 = fix_dihedrals_by_backbone_mapping(spc_1, spc_2, backbone_map=candidate) - fixed_spcs.append((fixed_spc_1, fixed_spc_2)) - backbone_1, backbone_2 = set(list(candidate.keys())), set(list(candidate.values())) - xyz1, xyz2 = fixed_spc_1.get_xyz(), fixed_spc_2.get_xyz() - xyz1 = xyz_from_data(coords=[xyz1['coords'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], - symbols=[xyz1['symbols'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], - isotopes=[xyz1['isotopes'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1]) - xyz2 = xyz_from_data(coords=[xyz2['coords'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], - symbols=[xyz2['symbols'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], - isotopes=[xyz2['isotopes'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2]) - no_gap_candidate = remove_gaps_from_values(candidate) - xyz2 = sort_xyz_using_indices(xyz2, indices=[v for k, v in sorted(no_gap_candidate.items(), - key=lambda item: item[0])]) - rmsds.append(compare_confs(xyz1=xyz1, xyz2=xyz2, rmsd_score=True)) - chosen_candidate_index = rmsds.index(min(rmsds)) - fixed_spc_1, fixed_spc_2 = fixed_spcs[chosen_candidate_index] - atom_map = map_hydrogens(fixed_spc_1, fixed_spc_2, candidate) - if map_type == 'list': - atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] - if atom_map is None and allow_backend_shift: - backend = 'QCElemental' - - if backend.lower() == 'qcelemental': - qcmol_1 = create_qc_mol(species=spc_1.copy()) - qcmol_2 = create_qc_mol(species=spc_2.copy()) - if qcmol_1 is None or qcmol_2 is None: + logger.warning(f'Could not identify superimposable candidates {spc_1} and {spc_2}.') return None - if len(qcmol_1.symbols) != len(qcmol_2.symbols): - raise ValueError(f'The number of atoms in spc1 ({spc_1.number_of_atoms}) must be equal ' - f'to the number of atoms in spc1 ({spc_2.number_of_atoms}).') - data = qcmol_2.align(ref_mol=qcmol_1, verbose=0, uno_cutoff=0.01) - atom_map = data[1]['mill'].atommap.tolist() - if map_type == 'dict': - atom_map = {key: val for key, val in enumerate(atom_map)} + if not len(candidates): + return None + rmsds, fixed_spcs = list(), list() + for candidate in candidates: + fixed_spc_1, fixed_spc_2 = fix_dihedrals_by_backbone_mapping(spc_1, spc_2, backbone_map=candidate) + fixed_spcs.append((fixed_spc_1, fixed_spc_2)) + backbone_1, backbone_2 = set(list(candidate.keys())), set(list(candidate.values())) + xyz1, xyz2 = fixed_spc_1.get_xyz(), fixed_spc_2.get_xyz() + xyz1 = xyz_from_data(coords=[xyz1['coords'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], + symbols=[xyz1['symbols'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], + isotopes=[xyz1['isotopes'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1]) + xyz2 = xyz_from_data(coords=[xyz2['coords'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], + symbols=[xyz2['symbols'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], + isotopes=[xyz2['isotopes'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2]) + no_gap_candidate = remove_gaps_from_values(candidate) + xyz2 = sort_xyz_using_indices(xyz2, indices=[v for k, v in sorted(no_gap_candidate.items(), key=lambda item: item[0])]) + rmsds.append(compare_confs(xyz1=xyz1, xyz2=xyz2, rmsd_score=True)) + chosen_candidate_index = rmsds.index(min(rmsds)) + fixed_spc_1, fixed_spc_2 = fixed_spcs[chosen_candidate_index] + atom_map = map_hydrogens(fixed_spc_1, fixed_spc_2, candidate) + if map_type == 'list': + atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] + return atom_map - if inc_vals is not None: - atom_map = [value + inc_vals for value in atom_map] + +def map_two_species_via_qcel(spc_1: Union[ARCSpecies, Species, Molecule], + spc_2: Union[ARCSpecies, Species, Molecule], + map_type: str = 'list', + ) -> Optional[Union[List[int], Dict[int, int]]]: + """ + Map the atoms in spc_1 to the atoms in spc_2 using the QCElemental method. + + Args: + spc_1 (Union[ARCSpecies, Species, Molecule]): Species 1. + spc_2 (Union[ARCSpecies, Species, Molecule]): Species 2. + map_type (str, optional): Whether to return a 'list' or a 'dict' map type. + + Returns: + Optional[Union[List[int], Dict[int, int]]]: + The atom maps. By default, a list of atom maps is returned. + """ + qcmol_1 = create_qc_mol(species=spc_1.copy()) + qcmol_2 = create_qc_mol(species=spc_2.copy()) + if qcmol_1 is None or qcmol_2 is None: + return None + if len(qcmol_1.symbols) != len(qcmol_2.symbols): + raise ValueError(f'The number of atoms in spc1 ({spc_1.number_of_atoms}) must be equal ' + f'to the number of atoms in spc1 ({spc_2.number_of_atoms}).') + data = qcmol_2.align(ref_mol=qcmol_1, verbose=0, uno_cutoff=0.01) + atom_map = data[1]['mill'].atommap.tolist() + if map_type == 'dict': + atom_map = {key: val for key, val in enumerate(atom_map)} return atom_map @@ -936,19 +976,20 @@ def check_atom_map(rxn) -> bool: Returns: bool Whether the atom mapping makes sense. """ - if len(rxn.atom_map) != sum([spc.number_of_atoms for spc in rxn.r_species]): - return False - r_elements, p_elements = list(), list() - for r_species in rxn.r_species: - r_elements.extend(list(r_species.get_xyz()['symbols'])) - for p_species in rxn.p_species: - p_elements.extend(list(p_species.get_xyz()['symbols'])) - for i, map_i in enumerate(rxn.atom_map): - if r_elements[i] != p_elements[map_i]: - break - for i in range(len(unique(rxn.atom_map))): - if unique(rxn.atom_map)[i] != i: + for atom_map in rxn.atom_maps: + if len(atom_map) != sum([spc.number_of_atoms for spc in rxn.r_species]): return False + r_elements, p_elements = list(), list() + for r_species in rxn.r_species: + r_elements.extend(list(r_species.get_xyz()['symbols'])) + for p_species in rxn.p_species: + p_elements.extend(list(p_species.get_xyz()['symbols'])) + for i, map_i in enumerate(atom_map): + if r_elements[i] != p_elements[map_i]: + break + for i in range(len(unique(atom_map))): + if unique(atom_map)[i] != i: + return False return True diff --git a/arc/reaction.py b/arc/reaction.py index 09e9dfbe42..29ac62b64a 100644 --- a/arc/reaction.py +++ b/arc/reaction.py @@ -82,9 +82,9 @@ class ARCReaction(object): ts_label (str): The :ref:`ARCSpecies ` label of the respective TS. preserve_param_in_scan (list): Entries are length two iterables of atom indices (1-indexed) between which distances and dihedrals of these pivots must be preserved. - atom_map (List[int]): An atom map, mapping the reactant atoms to the product atoms. - I.e., an atom map of [0, 2, 1] means that reactant atom 0 matches product atom 0, - reactant atom 1 matches product atom 2, and reactant atom 2 matches product atom 1. + atom_maps (List[List[int]]): Entries are atom maps, mapping the reactant atoms to the product atoms. + I.e., an atom map of [0, 2, 1] means that reactant atom 0 matches product atom 0, + reactant atom 1 matches product atom 2, and reactant atom 2 matches product atom 1. done_opt_r_n_p (bool): Whether the optimization of all reactants and products is complete. """ def __init__(self, @@ -117,7 +117,7 @@ def __init__(self, self.rmg_reactions = None self.ts_xyz_guess = ts_xyz_guess or xyz or list() self.preserve_param_in_scan = preserve_param_in_scan - self._atom_map = None + self._atom_maps = None self._charge = charge self._multiplicity = multiplicity if reaction_dict is not None: @@ -150,24 +150,24 @@ def __init__(self, self.check_atom_balance() @property - def atom_map(self): - """The reactants to products atom map""" - if self._atom_map is None \ + def atom_maps(self): + """The reactants to products atom maps""" + if self._atom_maps is None \ and all(species.get_xyz(generate=False) is not None for species in self.r_species + self.p_species): for backend in ["ARC", "QCElemental"]: - _atom_map = map_reaction(rxn=self, backend=backend) - if _atom_map is not None: - self._atom_map = _atom_map + _atom_maps = map_reaction(rxn=self, backend=backend) + if _atom_maps is not None: + self._atom_maps = _atom_maps break logger.error(f"The requested ARC reaction {self}, and it's reverse, could not be atom mapped using {backend}.") - if self._atom_map is None: + if self._atom_maps is None: logger.error(f"The requested ARC reaction {self} could not be atom mapped.") - return self._atom_map + return self._atom_maps - @atom_map.setter - def atom_map(self, value): + @atom_maps.setter + def atom_maps(self, value): """Allow setting the atom map""" - self._atom_map = value + self._atom_maps = value @property def charge(self): @@ -239,8 +239,8 @@ def as_dict(self, reaction_dict['p_species'] = [spc.as_dict(reset_atom_ids=reset_atom_ids) for spc in self.p_species] if self.ts_species is not None: reaction_dict['ts_species'] = self.ts_species.as_dict() - if self._atom_map is not None: - reaction_dict['atom_map'] = self._atom_map + if self._atom_maps is not None: + reaction_dict['atom_maps'] = self._atom_maps if self.done_opt_r_n_p is not None: reaction_dict['done_opt_r_n_p'] = self.done_opt_r_n_p if self.preserve_param_in_scan is not None: @@ -321,7 +321,7 @@ def from_dict(self, else reaction_dict['xyz'] if 'xyz' in reaction_dict else list() self.preserve_param_in_scan = reaction_dict['preserve_param_in_scan'] \ if 'preserve_param_in_scan' in reaction_dict else None - self.atom_map = reaction_dict['atom_map'] if 'atom_map' in reaction_dict else None + self.atom_maps = reaction_dict['atom_maps'] if 'atom_maps' in reaction_dict else None self.done_opt_r_n_p = reaction_dict['done_opt_r_n_p'] if 'done_opt_r_n_p' in reaction_dict else None def copy(self): @@ -342,7 +342,7 @@ def flip_reaction(self): ARCReaction: A copy of this object instance with flipped reactants and products. """ reaction_dict = self.as_dict(reset_atom_ids=True) - reset_keys = ['label', 'index', 'atom_map', 'rmg_reaction', + reset_keys = ['label', 'index', 'atom_maps', 'rmg_reaction', 'family', 'family_own_reverse', 'long_kinetic_description'] if 'r_species' in reaction_dict.keys() and 'p_species' in reaction_dict.keys(): reaction_dict['r_species'], reaction_dict['p_species'] = reaction_dict['p_species'], reaction_dict['r_species'] @@ -841,7 +841,9 @@ def get_number_of_atoms_in_reaction_zone(self) -> Optional[int]: labels = set(labels) return len(labels) - def get_single_mapped_product_xyz(self) -> Optional[ARCSpecies]: + def get_single_mapped_product_xyz(self, + am_index: int = 0, + ) -> Optional[ARCSpecies]: """ Get a copy of the product species with mapped cartesian coordinates of a reaction with a single product. @@ -852,7 +854,7 @@ def get_single_mapped_product_xyz(self) -> Optional[ARCSpecies]: logger.error(f'Can only return a mapped product for reactions with a single product, ' f'got {len(self.p_species)}.') return None - mapped_xyz = sort_xyz_using_indices(xyz_dict=self.p_species[0].get_xyz(), indices=self.atom_map) + mapped_xyz = sort_xyz_using_indices(xyz_dict=self.p_species[0].get_xyz(), indices=self.atom_maps[am_index]) mapped_product = ARCSpecies(label=self.p_species[0].label, mol=self.p_species[0].mol.copy(deep=True), multiplicity=self.p_species[0].multiplicity, @@ -921,7 +923,7 @@ def get_products_xyz(self, return_format='str') -> Union[dict, str]: xyz_dict['isotopes'] += xyz['isotopes'] xyz_dict['coords'] += xyz['coords'] xyz_dict = translate_to_center_of_mass(check_xyz_dict(xyz_dict)) - xyz_dict = sort_xyz_using_indices(xyz_dict=xyz_dict, indices=self.atom_map) + xyz_dict = sort_xyz_using_indices(xyz_dict=xyz_dict, indices=self.atom_maps[0]) if return_format == 'str': xyz_dict = xyz_to_str(xyz_dict) return xyz_dict diff --git a/arc/reaction_test.py b/arc/reaction_test.py index 93753ac62e..4d747c3aaa 100644 --- a/arc/reaction_test.py +++ b/arc/reaction_test.py @@ -374,7 +374,7 @@ def test_copy(self): self.assertEqual(self.rxn8.charge, rxn_copy.charge) self.assertTrue(check_atom_map(self.rxn8)) self.assertTrue(check_atom_map(rxn_copy)) - self.assertEqual(self.rxn8.atom_map, rxn_copy.atom_map) + self.assertEqual(self.rxn8.atom_maps, rxn_copy.atom_maps) self.assertEqual(tuple(spc.label for spc in self.rxn8.r_species), tuple(spc.label for spc in rxn_copy.r_species)) self.assertEqual(tuple(spc.label for spc in self.rxn8.p_species), @@ -403,10 +403,10 @@ def test_flip_reaction(self): self.assertEqual(self.rxn8.multiplicity, flipped_rxn.multiplicity) self.assertEqual(self.rxn8.charge, flipped_rxn.charge) self.assertTrue(check_atom_map(flipped_rxn)) - self.assertEqual(self.rxn8.atom_map[0], 0) - self.assertEqual(self.rxn8.atom_map[5], 4) - self.assertEqual(flipped_rxn.atom_map[0], 0) - self.assertEqual(flipped_rxn.atom_map[4], 5) + self.assertEqual(self.rxn8.atom_maps[0][0], 0) + self.assertEqual(self.rxn8.atom_maps[0][5], 4) + self.assertEqual(flipped_rxn.atom_maps[0][0], 0) + self.assertEqual(flipped_rxn.atom_maps[0][4], 5) self.assertEqual(tuple(spc.label for spc in self.rxn8.r_species), tuple(spc.label for spc in flipped_rxn.p_species)) self.assertEqual(tuple(spc.label for spc in self.rxn8.p_species), @@ -699,9 +699,8 @@ def test_get_number_of_atoms_in_reaction_zone(self): self.assertEqual(self.rxn10.get_number_of_atoms_in_reaction_zone(), 3) self.assertEqual(self.rxn11.get_number_of_atoms_in_reaction_zone(), 3) - def test_get_atom_map(self): - """Test getting an atom map for a reaction""" - + def test_get_atom_map_unimolecular(self): + """Test getting an atom map for unimolecular reactions""" # 1. trivial unimolecular: H2O <=> H2O h2o_xyz_1 = {'symbols': ('O', 'H', 'H'), 'isotopes': (16, 1, 1), 'coords': ((-0.0003283189391273643, 0.39781490416473486, 0.0), @@ -710,7 +709,7 @@ def test_get_atom_map(self): r_1 = ARCSpecies(label='H2O', smiles='O', xyz=h2o_xyz_1) p_1 = ARCSpecies(label='H2O', smiles='O', xyz=h2o_xyz_1) rxn_1 = ARCReaction(reactants=['H2O'], products=['H2O'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn_1.atom_map, [0, 1, 2]) + self.assertEqual(rxn_1.atom_maps[0], [0, 1, 2]) self.assertTrue(check_atom_map(rxn_1)) # 2. trivial unimolecular with an intentional mixed atom order: H2O <=> H2O @@ -720,10 +719,12 @@ def test_get_atom_map(self): (-0.0003283189391273643, 0.39781490416473486, 0.0))} p_1 = ARCSpecies(label='H2O', smiles='O', xyz=h2o_xyz_2) rxn_2 = ARCReaction(reactants=['H2O'], products=['H2O'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn_2.atom_map, [2, 0, 1]) + self.assertEqual(rxn_2.atom_maps[0], [2, 0, 1]) self.assertTrue(check_atom_map(rxn_2)) - # 3. trivial bimolecular: H + CH3NH2 <=> H2 + CH2NH2 + def test_get_atom_map_h_abstraction(self): + """Test getting an atom map for H abstraction reactions""" + # trivial: H + CH3NH2 <=> H2 + CH2NH2 ch3nh2_xyz = {'coords': ((-0.5734111454228507, 0.0203516083213337, 0.03088703933770556), (0.8105595891860601, 0.00017446498908627427, -0.4077728757313545), (-1.1234549667791063, -0.8123899006368857, -0.41607711106038836), @@ -749,21 +750,21 @@ def test_get_atom_map(self): r_species=[r_1, r_2], p_species=[p_1, p_2]) rxn_3.determine_family(self.rmgdb) self.assertEqual(rxn_3.family.label.lower(), "H_Abstraction".lower()) - self.assertIn(rxn_3.atom_map[0], [0, 1]) - self.assertEqual(rxn_3.atom_map[1:3], [2, 5]) + self.assertIn(rxn_3.atom_maps[0][0], [0, 1]) + self.assertEqual(rxn_3.atom_maps[0][1:3], [2, 5]) for index in [3, 4, 5]: - self.assertIn(rxn_3.atom_map[index], [0, 1, 3, 4]) - self.assertIn(rxn_3.atom_map[6:], [[7, 6], [6, 7]]) + self.assertIn(rxn_3.atom_maps[0][index], [0, 1, 3, 4]) + self.assertIn(rxn_3.atom_maps[0][6:], [[7, 6], [6, 7]]) self.assertTrue(check_atom_map(rxn_3)) - # 4. trivial bimolecular in reverse order: H + CH3NH2 <=> CH2NH2 + H2 + # trivial in reverse order: H + CH3NH2 <=> CH2NH2 + H2 rxn_4 = ARCReaction(reactants=['H', 'CH3NH2'], products=['CH2NH2', 'H2'], r_species=[r_1, r_2], p_species=[p_2, p_1]) - self.assertIn(rxn_4.atom_map[0], [6, 7]) - self.assertEqual(rxn_4.atom_map[1:3], [0, 3]) + self.assertIn(rxn_4.atom_maps[0][0], [6, 7]) + self.assertEqual(rxn_4.atom_maps[0][1:3], [0, 3]) for index in [3, 4, 5]: - self.assertIn(rxn_4.atom_map[index], [1, 2, 6, 7]) - self.assertIn(rxn_4.atom_map[6:], [[5, 4], [4, 5]]) + self.assertIn(rxn_4.atom_maps[0][index], [1, 2, 6, 7]) + self.assertIn(rxn_4.atom_maps[0][6:], [[5, 4], [4, 5]]) self.assertTrue(check_atom_map(rxn_4)) # H_Abstraction: C3H6O + C4H9O <=> C3H5O + C4H10O @@ -829,7 +830,7 @@ def test_get_atom_map(self): p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=c4h10o_xyz) rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] self.assertEqual(atom_map[0:4], [0, 1, 3, 4]) self.assertIn(atom_map[4], [5,6, 7]) self.assertIn(atom_map[5], [5,6, 7]) @@ -878,16 +879,18 @@ def test_get_atom_map(self): p_2 = ARCSpecies(label='N2H2(T)', smiles='[NH][NH]', xyz=n2h2_t_xyz) rxn = ARCReaction(reactants=['NH', 'N2H3'], products=['NH2', 'N2H2(T)'], r_species=[r_1, r_2], p_species=[p_1, p_2]) - self.assertEqual(rxn.atom_map[0], 0) - self.assertIn(rxn.atom_map[1], [1, 2]) - self.assertIn(rxn.atom_map[2], [3, 5]) - self.assertIn(rxn.atom_map[3], [3, 5]) - self.assertIn(rxn.atom_map[4], [4, 6]) - self.assertIn(rxn.atom_map[5], [1, 2, 4, 6]) - self.assertIn(rxn.atom_map[6], [1, 2, 4, 6]) - self.assertTrue(any(rxn.atom_map[index] in [1, 2] for index in [5, 6])) + self.assertEqual(rxn.atom_maps[0][0], 0) + self.assertIn(rxn.atom_maps[0][1], [1, 2]) + self.assertIn(rxn.atom_maps[0][2], [3, 5]) + self.assertIn(rxn.atom_maps[0][3], [3, 5]) + self.assertIn(rxn.atom_maps[0][4], [4, 6]) + self.assertIn(rxn.atom_maps[0][5], [1, 2, 4, 6]) + self.assertIn(rxn.atom_maps[0][6], [1, 2, 4, 6]) + self.assertTrue(any(rxn.atom_maps[0][index] in [1, 2] for index in [5, 6])) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_ring_opening(self): + """Test getting an atom map for ring opening reactions""" # Cyclopentadiene_scission: C6H6 <=> C6H6_2 c6h6_a_xyz = {'coords': ((1.465264096022479, 0.3555098886638667, 0.15268159347190322), (0.4583546746026421, 1.1352991023740606, -0.26555553330413073), @@ -932,9 +935,11 @@ def test_get_atom_map(self): 11 H u0 p0 c0 {4,S} 12 H u0 p0 c0 {5,S}""") rxn = ARCReaction(reactants=['C6H6_1'], products=['C6H6_b'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [1, 4, 2, 0, 5, 3, 10, 9, 8, 7, 6, 11]) + self.assertEqual(rxn.atom_maps[0], [1, 4, 2, 0, 5, 3, 10, 9, 8, 7, 6, 11]) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_disprop(self): + """Test getting an atom map for disproportionation reactions""" # Disproportionation: HO2 + NHOH <=> NH2OH + O2 nhoh_xyz = {'coords': ((0.5055094877826753, 0.03248552573561613, -0.443416250587286), (1.392367115364475, -0.021750569314658803, 0.07321920788090872), @@ -953,99 +958,13 @@ def test_get_atom_map(self): p_2 = ARCSpecies(label='NH2OH', smiles='NO', xyz=nh2oh_xyz) rxn = ARCReaction(reactants=['NHOH', 'HO2'], products=['O2', 'NH2OH'], r_species=[r_1, r_2], p_species=[p_1, p_2]) - self.assertEqual(rxn.atom_map[0], 2) + self.assertEqual(rxn.atom_maps[0][0], 2) for index in [1, 6]: - self.assertIn(rxn.atom_map[index], [4, 5]) - self.assertEqual(rxn.atom_map[2], 3) - self.assertEqual(rxn.atom_map[3], 6) + self.assertIn(rxn.atom_maps[0][index], [4, 5]) + self.assertEqual(rxn.atom_maps[0][2], 3) + self.assertEqual(rxn.atom_maps[0][3], 6) for index in [4, 5]: - self.assertIn(rxn.atom_map[index], [0, 1]) - self.assertTrue(check_atom_map(rxn)) - - # Intra_Disproportionation: C10H10_a <=> C10H10_b - c10h10_a_xyz = {'coords': ((3.1623638230700997, 0.39331289450005563, -0.031839117414963584), - (1.8784852381397288, 0.037685951926618944, -0.13659028131444134), - (0.9737380560194014, 0.5278617594060281, -1.1526858375270472), - (1.2607098516126556, 1.1809007875206383, -1.9621017164412065), - (-0.36396095305912823, -0.13214785064139675, -1.0200667625809143), - (-1.5172464644867296, 0.8364138939810618, -1.0669384323486588), - (-2.4922101649968655, 0.8316551483126366, -0.14124720277902958), - (-2.462598061982958, -0.09755474191953761, 0.9703503187569243), - (-1.4080417204047313, -0.8976377310686736, 1.1927020968566089), - (-0.27981087345916755, -0.8670643393461046, 0.29587765657632165), - (1.1395623815572733, -0.9147118621123697, 0.771368745020215), - (3.7901243915692864, -0.006544237180536178, 0.7580206603561134), - (3.6186251824572455, 1.0920401631166292, -0.725695658374561), - (-0.4799044636709365, -0.8577283498506146, -1.8345168113636874), - (-1.5704890060131314, 1.527002009812866, -1.902575985299536), - (-3.3260277144990296, 1.5238536460491903, -0.20338465526703625), - (-3.311126364299293, -0.09969554359088921, 1.6478137927333953), - (-1.3707042898204835, -1.549541647625315, 2.0589774409040964), - (1.5338362221707007, -1.9310023570889727, 0.6663504223502944), - (1.2246749300961473, -0.5970975942012858, 1.816181327157103)), - 'isotopes': (12, 12, 12, 1, 12, 12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1, 1, 1, 1), - 'symbols': ('C', 'C', 'C', 'H', 'C', 'C', 'C', 'C', 'C', 'C', 'C', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H')} - c10h10_b_xyz = {'coords': ((3.247237794328524, -0.13719671162966918, 0.19555833918937052), - (1.9094861712282774, -0.08067655688828143, 0.14898941432495702), - (0.9973729357914858, -1.2703386896415134, -0.09322848415119056), - (-0.37904449715218924, -0.6747166782032148, -0.049044448345326556), - (-0.32906812544096026, 0.704634441388649, 0.189424753183012), - (-1.4900181263846768, 1.4572613706024167, 0.2695747550348709), - (-2.715200996994148, 0.8069241052920498, 0.10660938013945513), - (-2.765284083663716, -0.5753713833636181, -0.13236922431004927), - (-1.5909002849280705, -1.3270914347507115, -0.21179882275795825), - (1.0862366144301145, 1.1823049698313937, 0.33079658088902575), - (3.8424769924852367, 0.7530758608805569, 0.37314678191170336), - (3.7762437608797406, -1.0749685445597326, 0.05710603017340202), - (1.1128196175313243, -2.0170485762246773, 0.6986324476157837), - (1.187449599052061, -1.7129398667445945, -1.0760419644685346), - (-1.453108430051206, 2.525963604437891, 0.45426129138400156), - (-3.639988653002051, 1.3756767310587803, 0.16518163487425436), - (-3.7283956370857467, -1.0643593255501977, -0.2566648708585298), - (-1.631427244782937, -2.3956407728893367, -0.3966116183664473), - (1.3188711462571718, 1.9143096670969255, -0.4489453399950017), - (1.2442414475018486, 1.6101977898569013, 1.3257284397785851)), - 'isotopes': (12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), - 'symbols': ('C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H')} - r_1 = ARCSpecies(label='C10H10_a', smiles='C=C1[CH]C2C=CC=C[C]2C1', xyz=c10h10_a_xyz, multiplicity=3) - p_1 = ARCSpecies(label='C10H10_b', smiles='C=C1CC2=C(C=CC=C2)C1', xyz=c10h10_b_xyz) - rxn = ARCReaction(reactants=['C10H10_a'], products=['C10H10_b'], - r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [0, 1, 8, 13, 3, 2, 7, 6, 5, 4, 9, 10, 17, 12, 11, 16, 15, 14, 19, 18]) - self.assertTrue(check_atom_map(rxn)) - - # intra_NO2_ONO_conversion: C2H5NO2 <=> C2H5ONO - c6h5_xyz = {'coords': ((1.8953828083622057, 0.8695975650550358, 0.6461465212661076), - (1.3601473931706598, -0.04212583715410005, 0.0034200061443233247), - (1.8529583069008781, -0.6310931351538215, -0.9666668585141432), - (-0.010154355673379136, -0.4652844276756663, 0.43320585211058743), - (-1.0281604639422022, 0.36855062612122236, -0.3158851121891869), - (-0.11071296591935365, -1.5314728469286516, 0.20909234121344752), - (-0.07635985361458197, -0.31625218083177237, 1.5151037167736001), - (-2.042322710601489, 0.08102183703582924, -0.021667016484293297), - (-0.9033569412063314, 1.436005790671757, -0.10388682333330314), - (-0.937421217476434, 0.23105260886017234, -1.3988626269871478)), - 'isotopes': (16, 14, 16, 12, 12, 1, 1, 1, 1, 1), - 'symbols': ('O', 'N', 'O', 'C', 'C', 'H', 'H', 'H', 'H', 'H')} - c7h5o_xyz = {'coords': ((-1.3334725178745668, 0.2849178019354427, 0.4149005134933577), - (-0.08765353373275289, 0.24941420749682627, -0.4497882845360618), - (1.0488580188184402, 0.3986394744609146, 0.39515448276833964), - (2.2292240798482883, 0.36629637181188207, -0.4124684043339001), - (3.2413605054484185, 0.4928521621538312, 0.283008378837631), - (-1.3088339518827734, -0.5173661350567303, 1.1597967522753032), - (-2.23462275856269, 0.17332354052924734, -0.19455307765792382), - (-1.393828440234405, 1.2294860794610234, 0.9656140588162426), - (-0.12370667081323389, 1.0672740524773998, -1.1795070012935482), - (-0.037324731014725374, -0.7080479312151163, -0.9821574183694773)), - 'isotopes': (12, 12, 16, 14, 16, 1, 1, 1, 1, 1), - 'symbols': ('C', 'C', 'O', 'N', 'O', 'H', 'H', 'H', 'H', 'H')} - r_1 = ARCSpecies(label='C2H5NO2', smiles='[O-][N+](=O)CC', xyz=c6h5_xyz) - p_1 = ARCSpecies(label='C2H5ONO', smiles='CCON=O', xyz=c7h5o_xyz) - rxn = ARCReaction(reactants=['C2H5NO2'], products=['C2H5ONO'], - r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [4, 3, 2, 1, 0, 8, 9, 6, 5, 7]) + self.assertIn(rxn.atom_maps[0][index], [0, 1]) self.assertTrue(check_atom_map(rxn)) # Disproportionation: C4H7 + O2 <=> HO2 + C4H6 @@ -1080,11 +999,11 @@ def test_get_atom_map(self): p_2 = ARCSpecies(label='C4H6', smiles='C=CC=C', xyz=c4h6_xyz) rxn = ARCReaction(reactants=['C4H7', 'O2'], products=['HO2', 'C4H6'], r_species=[r_1, r_2], p_species=[p_1, p_2]) - self.assertEqual(rxn.atom_map[:5], [3, 4, 5, 10, 6]) - self.assertIn(rxn.atom_map[5:7], [[8, 7], [7, 8]]) - self.assertEqual(rxn.atom_map[7], 9) - self.assertIn(tuple(rxn.atom_map[8:11]), permutations([2, 11, 12])) - self.assertIn(tuple(rxn.atom_map[11:]), permutations([0, 1])) + self.assertEqual(rxn.atom_maps[0][:5], [3, 4, 5, 10, 6]) + self.assertIn(rxn.atom_maps[0][5:7], [[8, 7], [7, 8]]) + self.assertEqual(rxn.atom_maps[0][7], 9) + self.assertIn(tuple(rxn.atom_maps[0][8:11]), permutations([2, 11, 12])) + self.assertIn(tuple(rxn.atom_maps[0][11:]), permutations([0, 1])) # Disproportionation: HO2 + NHOH <=> NH2OH + O2 nhoh_xyz = {'coords': ((0.5055094877826753, 0.03248552573561613, -0.443416250587286), @@ -1104,13 +1023,13 @@ def test_get_atom_map(self): p_2 = ARCSpecies(label='NH2OH', smiles='NO', xyz=nh2oh_xyz) rxn = ARCReaction(reactants=['NHOH', 'HO2'], products=['O2', 'NH2OH'], r_species=[r_1, r_2], p_species=[p_1, p_2]) - self.assertEqual(rxn.atom_map[0], 2) + self.assertEqual(rxn.atom_maps[0][0], 2) for index in [1, 6]: - self.assertIn(rxn.atom_map[index], [4, 5]) - self.assertEqual(rxn.atom_map[2], 3) - self.assertEqual(rxn.atom_map[3], 6) + self.assertIn(rxn.atom_maps[0][index], [4, 5]) + self.assertEqual(rxn.atom_maps[0][2], 3) + self.assertEqual(rxn.atom_maps[0][3], 6) for index in [4, 5]: - self.assertIn(rxn.atom_map[index], [0, 1]) + self.assertIn(rxn.atom_maps[0][index], [0, 1]) self.assertTrue(check_atom_map(rxn)) # Intra_Disproportionation: C10H10_a <=> C10H10_b @@ -1164,57 +1083,66 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C10H10_b', smiles='C=C1CC2=C(C=CC=C2)C1', xyz=c10h10_b_xyz) rxn = ARCReaction(reactants=['C10H10_a'], products=['C10H10_b'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [0, 1, 8, 13, 3, 2, 7, 6, 5, 4, 9, 10, 17, 12, 11, 16, 15, 14, 19, 18]) + self.assertEqual(rxn.atom_maps[0], [0, 1, 8, 13, 3, 2, 7, 6, 5, 4, 9, 10, 17, 12, 11, 16, 15, 14, 19, 18]) self.assertTrue(check_atom_map(rxn)) - # Cyclopentadiene_scission: C6H6 <=> C6H6_2 - c6h6_a_xyz = {'coords': ((1.465264096022479, 0.3555098886638667, 0.15268159347190322), - (0.4583546746026421, 1.1352991023740606, -0.26555553330413073), - (-0.7550043760214846, 0.35970165318809594, -0.5698935045151712), - (-1.485327813119871, -0.35660657095915016, 0.46119177830578917), - (-0.3414477960946828, -1.060779229397218, -0.11686056681841692), - (0.9879417277856641, -1.006839916409751, 0.12489717473407935), - (2.4630837864551887, 0.6629994259328668, 0.4197578798464181), - (0.5110882588097015, 2.2100951208919897, -0.3734820378556644), - (-1.1192886361027838, 0.384286081689225, -1.5897813181530946), - (-2.453224961870327, -0.7758708758357847, 0.2158838729688473), - (-1.3859013659398718, -0.054382091828296085, 1.4971154213962072), - (1.6544624054733257, -1.8534125883098933, 0.0440452399232336)), - 'isotopes': (12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1), - 'symbols': ('C', 'C', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H')} - c6h6_b_xyz = {'coords': ((-1.474267041853848, 0.27665693719971857, -0.31815898666696507), - (-0.25527025747758825, 1.1936776717612125, -0.2432148642540069), - (0.9917471212521393, 0.7578589393970138, 0.059037260524552534), - (1.2911962562420976, -0.6524103892231805, 0.34598643264742923), - (0.321535921890914, -1.5867102018006056, 0.32000545365633654), - (-0.9417407846918554, -1.043897260224426, -0.002820356559266387), - (-2.2262364004658077, 0.5956762298613206, 0.40890113659975075), - (-1.90597332290244, 0.31143075666839354, -1.3222845692785703), - (-0.4221153027089989, 2.2469871640348815, -0.4470234892644997), - (1.824518548011024, 1.4543788790156666, 0.0987362566117616), - (2.3174577767359237, -0.9162726684959432, 0.5791638390925197), - (0.4791474859684761, -2.637376058194065, 0.5216718868909702)), - 'isotopes': (12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1), - 'symbols': ('C', 'C', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H')} - r_1 = ARCSpecies(label='C6H6_1', smiles='C1=CC2CC2=C1', xyz=c6h6_a_xyz) - p_1 = ARCSpecies(label='C6H6_b', xyz=c6h6_b_xyz, adjlist="""multiplicity 1 - 1 C u0 p0 c0 {2,S} {6,S} {7,S} {8,S} - 2 C u0 p0 c0 {1,S} {3,D} {9,S} - 3 C u0 p0 c0 {2,D} {4,S} {10,S} - 4 C u0 p0 c0 {3,S} {5,D} {11,S} - 5 C u0 p0 c0 {4,D} {6,S} {12,S} - 6 C u0 p1 c0 {1,S} {5,S} - 7 H u0 p0 c0 {1,S} - 8 H u0 p0 c0 {1,S} - 9 H u0 p0 c0 {2,S} - 10 H u0 p0 c0 {3,S} - 11 H u0 p0 c0 {4,S} - 12 H u0 p0 c0 {5,S}""") - rxn = ARCReaction(reactants=['C6H6_1'], products=['C6H6_b'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [1, 4, 2, 0, 5, 3, 10, 9, 8, 7, 6, 11]) + # Intra_Disproportionation: C10H10_a <=> C10H10_b + c10h10_a_xyz = {'coords': ((3.1623638230700997, 0.39331289450005563, -0.031839117414963584), + (1.8784852381397288, 0.037685951926618944, -0.13659028131444134), + (0.9737380560194014, 0.5278617594060281, -1.1526858375270472), + (1.2607098516126556, 1.1809007875206383, -1.9621017164412065), + (-0.36396095305912823, -0.13214785064139675, -1.0200667625809143), + (-1.5172464644867296, 0.8364138939810618, -1.0669384323486588), + (-2.4922101649968655, 0.8316551483126366, -0.14124720277902958), + (-2.462598061982958, -0.09755474191953761, 0.9703503187569243), + (-1.4080417204047313, -0.8976377310686736, 1.1927020968566089), + (-0.27981087345916755, -0.8670643393461046, 0.29587765657632165), + (1.1395623815572733, -0.9147118621123697, 0.771368745020215), + (3.7901243915692864, -0.006544237180536178, 0.7580206603561134), + (3.6186251824572455, 1.0920401631166292, -0.725695658374561), + (-0.4799044636709365, -0.8577283498506146, -1.8345168113636874), + (-1.5704890060131314, 1.527002009812866, -1.902575985299536), + (-3.3260277144990296, 1.5238536460491903, -0.20338465526703625), + (-3.311126364299293, -0.09969554359088921, 1.6478137927333953), + (-1.3707042898204835, -1.549541647625315, 2.0589774409040964), + (1.5338362221707007, -1.9310023570889727, 0.6663504223502944), + (1.2246749300961473, -0.5970975942012858, 1.816181327157103)), + 'isotopes': (12, 12, 12, 1, 12, 12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1, 1, 1, 1), + 'symbols': ('C', 'C', 'C', 'H', 'C', 'C', 'C', 'C', 'C', 'C', 'C', + 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H')} + c10h10_b_xyz = {'coords': ((3.247237794328524, -0.13719671162966918, 0.19555833918937052), + (1.9094861712282774, -0.08067655688828143, 0.14898941432495702), + (0.9973729357914858, -1.2703386896415134, -0.09322848415119056), + (-0.37904449715218924, -0.6747166782032148, -0.049044448345326556), + (-0.32906812544096026, 0.704634441388649, 0.189424753183012), + (-1.4900181263846768, 1.4572613706024167, 0.2695747550348709), + (-2.715200996994148, 0.8069241052920498, 0.10660938013945513), + (-2.765284083663716, -0.5753713833636181, -0.13236922431004927), + (-1.5909002849280705, -1.3270914347507115, -0.21179882275795825), + (1.0862366144301145, 1.1823049698313937, 0.33079658088902575), + (3.8424769924852367, 0.7530758608805569, 0.37314678191170336), + (3.7762437608797406, -1.0749685445597326, 0.05710603017340202), + (1.1128196175313243, -2.0170485762246773, 0.6986324476157837), + (1.187449599052061, -1.7129398667445945, -1.0760419644685346), + (-1.453108430051206, 2.525963604437891, 0.45426129138400156), + (-3.639988653002051, 1.3756767310587803, 0.16518163487425436), + (-3.7283956370857467, -1.0643593255501977, -0.2566648708585298), + (-1.631427244782937, -2.3956407728893367, -0.3966116183664473), + (1.3188711462571718, 1.9143096670969255, -0.4489453399950017), + (1.2442414475018486, 1.6101977898569013, 1.3257284397785851)), + 'isotopes': (12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), + 'symbols': ('C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', + 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H')} + r_1 = ARCSpecies(label='C10H10_a', smiles='C=C1[CH]C2C=CC=C[C]2C1', xyz=c10h10_a_xyz, multiplicity=3) + p_1 = ARCSpecies(label='C10H10_b', smiles='C=C1CC2=C(C=CC=C2)C1', xyz=c10h10_b_xyz) + rxn = ARCReaction(reactants=['C10H10_a'], products=['C10H10_b'], + r_species=[r_1], p_species=[p_1]) + self.assertEqual(rxn.atom_maps[0], [0, 1, 8, 13, 3, 2, 7, 6, 5, 4, 9, 10, 17, 12, 11, 16, 15, 14, 19, 18]) self.assertTrue(check_atom_map(rxn)) - # intra_NO2_ONO_conversion: C2H5NO2 <=> C2H5ONO + def test_get_atom_map_intra_no2_ono_conversion(self): + """Test getting an atom map for intra_NO2_ONO_conversion reactions""" + # C2H5NO2 <=> C2H5ONO c6h5_xyz = {'coords': ((1.8953828083622057, 0.8695975650550358, 0.6461465212661076), (1.3601473931706598, -0.04212583715410005, 0.0034200061443233247), (1.8529583069008781, -0.6310931351538215, -0.9666668585141432), @@ -1243,11 +1171,12 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C2H5ONO', smiles='CCON=O', xyz=c7h5o_xyz) rxn = ARCReaction(reactants=['C2H5NO2'], products=['C2H5ONO'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [4, 3, 2, 1, 0, 8, 9, 6, 5, 7]) + self.assertEqual(rxn.atom_maps[0], [4, 3, 2, 1, 0, 8, 9, 6, 5, 7]) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_birad_to_alkene(self): + """Test getting an atom map for birad to alkene reactions""" # 1,2-Birad_to_alkene: SO2(T) => SO2(S) - # This test fails due to a problem in xyz perception of SO2(T). so2_t_xyz = {'coords': ((0.02724478716956233, 0.6093829407458188, 0.0), (-1.3946381818031768, -0.24294788636871906, 0.0), (1.3673933946336125, -0.36643505437710233, 0.0)), @@ -1260,8 +1189,11 @@ def test_get_atom_map(self): r_1 = ARCSpecies(label='SO2(T)', smiles='O=[S][O]', multiplicity=3, xyz=so2_t_xyz) p_1 = ARCSpecies(label='SO2(S)', smiles='O=S=O', multiplicity=1, xyz=so2_s_xyz) rxn = ARCReaction(reactants=['SO2(T)'], products=['SO2(S)'], r_species=[r_1], p_species=[p_1]) - self.assertEqual(rxn.atom_map[0], 1) + self.assertEqual(rxn.atom_maps[0][0], 1) self.assertTrue(check_atom_map(rxn)) + + def test_get_atom_map_birad_r_rec(self): + """Test getting an atom map for birad R recombination reactions""" # F[C]F + [CH3] <=> F[C](F)C r1_xyz = {'symbols': ('F', 'C', 'F'), 'isotopes': (19, 12, 19), @@ -1286,12 +1218,14 @@ def test_get_atom_map(self): rxn = ARCReaction(r_species=[ARCSpecies(label="r1", smiles="F[C]F", xyz=r1_xyz), ARCSpecies(label="r2", smiles="[CH3]", xyz=r2_xyz)], p_species=[ARCSpecies(label="p1", smiles="F[C](F)C", xyz=p1_xyz)]) - self.assertIn(rxn.atom_map[:2], [[0, 1], [1, 0]]) - self.assertEqual(rxn.atom_map[2], 2) - self.assertEqual(rxn.atom_map[3], 3) - self.assertIn(tuple(rxn.atom_map[4:]), list(permutations([4, 5, 6]))) + self.assertIn(rxn.atom_maps[0][:2], [[0, 1], [1, 0]]) + self.assertEqual(rxn.atom_maps[0][2], 2) + self.assertEqual(rxn.atom_maps[0][3], 3) + self.assertIn(tuple(rxn.atom_maps[0][4:]), list(permutations([4, 5, 6]))) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_nh3_elimination(self): + """Test getting an atom map for NH3 elimination reactions""" # 1,2_NH3_elimination: NCC <=> C2H4 + NH3 ncc_xyz = {'coords': ((1.1517341397735719, -0.37601689454792764, -0.5230788502681245), (0.2893395715754821, 0.449973844025586, 0.3114935868175311), @@ -1324,12 +1258,14 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C2H4', smiles='C=C', xyz=c2h4_xyz) p_2 = ARCSpecies(label='NH3', smiles='N', xyz=nh3_xyz) rxn = ARCReaction(reactants=['NCC'], products=['C2H4', 'NH3'], r_species=[r_1], p_species=[p_1, p_2]) - self.assertEqual(rxn.atom_map[0], 6) - self.assertIn(rxn.atom_map[1:3], [[1, 0], [0, 1]]) - self.assertIn(tuple(rxn.atom_map[3:5]+[rxn.atom_map[7]]), permutations([7, 8, 9])) - self.assertIn(tuple(rxn.atom_map[5:7]+rxn.atom_map[8:]), permutations([2, 3, 4, 5])) + self.assertEqual(rxn.atom_maps[0][0], 6) + self.assertIn(rxn.atom_maps[0][1:3], [[1, 0], [0, 1]]) + self.assertIn(tuple(rxn.atom_maps[0][3:5]+[rxn.atom_maps[0][7]]), permutations([7, 8, 9])) + self.assertIn(tuple(rxn.atom_maps[0][5:7]+rxn.atom_maps[0][8:]), permutations([2, 3, 4, 5])) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_cycloaddition(self): + """Test getting an atom map for Cycloaddition reactions""" # 1+2_Cycloaddition: CH2 + C2H4 <=> C3H6 c2h4_xyz = {'coords': ((0.6664040429179742, 0.044298334171779405, -0.0050238049104911735), (-0.6664040438461246, -0.04429833352898575, 0.00502380522486473), @@ -1356,11 +1292,13 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='cC3H6', smiles='C1CC1', xyz=c_c3h6_xyz) rxn = ARCReaction(reactants=['CH2', 'C2H4'], products=['cC3H6'], r_species=[r_1, r_2], p_species=[p_1]) for index in [0, 3, 4]: - self.assertIn(rxn.atom_map[index], [0, 1, 2]) + self.assertIn(rxn.atom_maps[0][index], [0, 1, 2]) for index in [1, 2, 5, 6, 7, 8]: - self.assertIn(rxn.atom_map[index], [3, 4, 5, 6, 7, 8]) + self.assertIn(rxn.atom_maps[0][index], [3, 4, 5, 6, 7, 8]) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_insertion(self): + """Test getting an atom map for insertion reactions""" # 1,2_Insertion_CO: C4H10 + CO <=> C5H10O c4h10_xyz = {'coords': ((-0.5828455298013108, 1.3281531294599287, -0.04960015063595639), (0.20452033859928953, 0.05503751610159247, -0.351590668388836), @@ -1401,17 +1339,17 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C5H10O', smiles='CC(C)(C)C=O', xyz=c5h10o_xyz) rxn = ARCReaction(reactants=['C4H10', 'CO'], products=['C5H10O'], r_species=[r_1, r_2], p_species=[p_1]) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] self.assertEqual(atom_map[:4], [0, 1, 2, 3]) - self.assertIn(tuple(rxn.atom_map[4:7]), permutations([6, 7, 8])) + self.assertIn(tuple(atom_map[4:7]), permutations([6, 7, 8])) self.assertEqual(atom_map[7], 15) - self.assertIn(tuple(rxn.atom_map[8:11]), permutations([9, 10, 11])) - self.assertIn(tuple(rxn.atom_map[11:14]), permutations([12, 13, 14])) + self.assertIn(tuple(atom_map[8:11]), permutations([9, 10, 11])) + self.assertIn(tuple(atom_map[11:14]), permutations([12, 13, 14])) self.assertEqual(atom_map[14:], [4, 5]) self.assertTrue(check_atom_map(rxn)) # same reaction in reverse: rxn_rev = ARCReaction(r_species=[p_1], p_species=[r_1, r_2]) - atom_map = rxn_rev.atom_map + atom_map = rxn_rev.atom_maps[0] for index in [0, 2, 3]: self.assertIn(atom_map[index], [0, 2, 3]) self.assertEqual(atom_map[1], 1) @@ -1420,6 +1358,8 @@ def test_get_atom_map(self): self.assertEqual(atom_map[15], 7) self.assertTrue(check_atom_map(rxn_rev)) + def test_get_atom_map_diels_alder(self): + """Test getting an atom map for Diels Alder reactions""" # Diels_alder_addition: C5H8 + C6H10 <=> C11H18 c5h8_xyz = {'coords': ((2.388426506127341, -0.6020682478448856, -0.8986239521455471), (1.396815470095451, 0.2559764141247285, -0.632876393172657), @@ -1491,18 +1431,20 @@ def test_get_atom_map(self): r_2 = ARCSpecies(label='C6H10', smiles='CC=CC=CC', xyz=c6h10_xyz) p_1 = ARCSpecies(label='C11H18', smiles='C=CC1C(C)C=CC(C)C1C', xyz=c11h18_xyz) rxn = ARCReaction(reactants=['C5H8', 'C6H10'], products=['C11H18'], r_species=[r_1, r_2], p_species=[p_1]) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] self.assertEqual(atom_map[:5], [0, 1, 2, 9, 10]) self.assertEqual(atom_map[:5], [0, 1, 2, 9, 10]) self.assertIn(atom_map[5:7], [[11, 12], [12, 11]]) self.assertEqual(atom_map[7:10], [13, 14, 25]) - self.assertIn(tuple(rxn.atom_map[10:13]), permutations([26, 27, 28])) + self.assertIn(tuple(atom_map[10:13]), permutations([26, 27, 28])) self.assertEqual(atom_map[13:19], [8, 7, 6, 5, 3, 4]) - self.assertIn(tuple(rxn.atom_map[19:22]), permutations([24, 23, 22])) + self.assertIn(tuple(atom_map[19:22]), permutations([24, 23, 22])) self.assertEqual(atom_map[22:26], [21, 20, 19, 15]) - self.assertIn(tuple(rxn.atom_map[26:]), permutations([16, 17, 18])) + self.assertIn(tuple(atom_map[26:]), permutations([16, 17, 18])) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_add_endocyclic(self): + """Test getting an atom map for add endocyclic reactions""" # Intra_R_Add_Endocyclic: C9H15_a <=> C9H15_b c9h15_a_xyz = {'coords': ((3.2994642637411093, -0.9763218631003405, -0.6681519125224107), (2.092867397835492, -0.585345209944081, -1.094234941414971), @@ -1562,17 +1504,19 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C9H15_b', smiles='C=CC1(C)C[C](C)C1C', xyz=c9h15_b_xyz) rxn = ARCReaction(reactants=['C9H15_a'], products=['C9H15_b'], r_species=[r_1], p_species=[p_1]) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] self.assertEqual(atom_map[0:9], list(range(9))) self.assertIn(atom_map[9:11], [[9, 10], [10, 9]]) self.assertEqual(atom_map[11], 11) - self.assertIn(tuple(rxn.atom_map[12:15]), permutations([13, 14, 12])) + self.assertIn(tuple(atom_map[12:15]), permutations([13, 14, 12])) self.assertIn(atom_map[15:17], [[15, 16], [16, 15]]) - self.assertIn(tuple(rxn.atom_map[17:20]), permutations([18, 17, 19])) + self.assertIn(tuple(atom_map[17:20]), permutations([18, 17, 19])) self.assertEqual(atom_map[20], 20) - self.assertIn(tuple(rxn.atom_map[21:]), permutations([23, 21, 22])) + self.assertIn(tuple(atom_map[21:]), permutations([23, 21, 22])) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_addition(self): + """Test getting an atom map for addition reactions""" # R_Addition_COm: C6H5 + CO <=> C7H5O c6h5_xyz = {'coords': ((0.1817676212163122, -1.6072341699404684, -0.014610043584505971), (1.3027386938520413, -0.7802649159703986, -0.0076490984025043415), @@ -1607,20 +1551,22 @@ def test_get_atom_map(self): p_1 = ARCSpecies(label='C7H5O', smiles='O=[C]c1ccccc1', xyz=c7h5o_xyz) rxn = ARCReaction(reactants=['C6H5', 'CO'], products=['C7H5O'], r_species=[r_1, r_2], p_species=[p_1]) - self.assertEqual(rxn.atom_map, [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 0]) + self.assertEqual(rxn.atom_maps[0], [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 0]) self.assertTrue(check_atom_map(rxn)) + def test_get_atom_map_keto_enol(self): + """Test getting an atom map for keto-enol reactions""" # Keto-enol isomerization: ch2choh <=> ch3cho r_1 = ARCSpecies(label='CH2CHOH', smiles='C=CO', xyz=self.ch2choh_xyz) p_1 = ARCSpecies(label='CH3CHO', smiles='CC=O', xyz=self.ch3cho_xyz) rxn = ARCReaction(r_species=[r_1], p_species=[p_1]) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] self.assertTrue(check_atom_map(rxn)) self.assertTrue(atom_map[:3], [0, 1, 2]) self.assertIn(tuple(atom_map[3:5]+[atom_map[-1]]), permutations([3, 4, 5])) self.assertEqual(atom_map[5], 6) - # Same species in products + # Same species in products rxn = ARCReaction(r_species=[ARCSpecies(label="r", smiles = 'C=C[CH]CC[CH]C=C')], p_species=[ARCSpecies(label="p1", smiles= 'C=CC=C'), ARCSpecies(label="p2", smiles= 'C=CC=C')]) @@ -1631,7 +1577,7 @@ def test_get_atom_map(self): h_symmetry2 = [4, 5, 8, 9, 14 ,15, 18, 19] # symmetry of hydrogens with two other hydrogen second order neighbor label_species_atoms(rxn.r_species) - atom_map = rxn.atom_map + atom_map = rxn.atom_maps[0] for atom in rxn.r_species[0].mol.copy(deep=True).atoms: if atom.symbol == "C": relevant_atom = atom @@ -1743,11 +1689,11 @@ def test_get_reactants_xyz(self): p_1 = ARCSpecies(label='C2H4O', smiles='CC=O', xyz=c2h4o_xyz) p_2 = ARCSpecies(label='HO2', smiles='O[O]', xyz=ho2_xyz) rxn = ARCReaction(r_species=[r_1], p_species=[p_1, p_2]) - self.assertIn(rxn.atom_map[0:5], [[0, 1, 2, 8, 7], [0, 1, 2, 7, 8]]) + self.assertIn(rxn.atom_maps[0][0:5], [[0, 1, 2, 8, 7], [0, 1, 2, 7, 8]]) for index in [5, 6, 7]: - self.assertIn(rxn.atom_map[index], [3, 4, 5]) - self.assertEqual(rxn.atom_map[8], 6) - self.assertEqual(rxn.atom_map[9], 9) + self.assertIn(rxn.atom_maps[0][index], [3, 4, 5]) + self.assertEqual(rxn.atom_maps[0][8], 6) + self.assertEqual(rxn.atom_maps[0][9], 9) self.assertTrue(check_atom_map(rxn)) def test_get_single_mapped_product_xyz(self): @@ -1764,7 +1710,7 @@ def test_get_single_mapped_product_xyz(self): rxn_1 = ARCReaction(reactants=['H2O'], products=['H2O'], r_species=[r_1], p_species=[p_1]) mapped_product = rxn_1.get_single_mapped_product_xyz() - self.assertEqual(rxn_1.atom_map, [2, 0, 1]) + self.assertEqual(rxn_1.atom_maps[0], [2, 0, 1]) self.assertTrue(check_atom_map(rxn_1)) expected_xyz = {'symbols': ('O', 'H', 'H'), 'isotopes': (16, 1, 1), 'coords': ((-0.00032832, 0.3978149, 0.0), (-0.76330345, -0.19953755, 0.0), @@ -1794,10 +1740,10 @@ def test_get_single_mapped_product_xyz(self): rxn_2 = ARCReaction(r_species=[reactant], p_species=[product]) self.assertTrue(check_atom_map(rxn_2)) mapped_product = rxn_2.get_single_mapped_product_xyz() - self.assertEqual(rxn_2.atom_map[:6], [0, 1, 2, 3, 4, 5]) - self.assertIn(rxn_2.atom_map[6], [6, 8]) - self.assertIn(rxn_2.atom_map[7], [6, 7]) - self.assertIn(rxn_2.atom_map[8], [7, 8]) + self.assertEqual(rxn_2.atom_maps[0][:6], [0, 1, 2, 3, 4, 5]) + self.assertIn(rxn_2.atom_maps[0][6], [6, 8]) + self.assertIn(rxn_2.atom_maps[0][7], [6, 7]) + self.assertIn(rxn_2.atom_maps[0][8], [7, 8]) expected_xyz = {'symbols': ('C', 'C', 'N', 'O', 'N', 'N', 'H', 'H', 'H'), 'isotopes': (12, 12, 14, 16, 14, 14, 1, 1, 1), 'coords': ((-1.0108, -0.0114, -0.061), (0.478, 0.0191, 0.0139), (1.2974, -0.993, 0.4693), @@ -1828,10 +1774,10 @@ def test_get_single_mapped_product_xyz(self): rxn_2 = ARCReaction(r_species=[reactant], p_species=[product]) self.assertTrue(check_atom_map(rxn_2)) mapped_product = rxn_2.get_single_mapped_product_xyz() - self.assertEqual(rxn_2.atom_map[:6], [0, 1, 2, 3, 4, 5]) - self.assertIn(rxn_2.atom_map[6], [6, 8]) - self.assertIn(rxn_2.atom_map[7], [6, 7]) - self.assertIn(rxn_2.atom_map[8], [7, 8]) + self.assertEqual(rxn_2.atom_maps[0][:6], [0, 1, 2, 3, 4, 5]) + self.assertIn(rxn_2.atom_maps[0][6], [6, 8]) + self.assertIn(rxn_2.atom_maps[0][7], [6, 7]) + self.assertIn(rxn_2.atom_maps[0][8], [7, 8]) expected_xyz = {'symbols': ('C', 'C', 'N', 'O', 'N', 'N', 'H', 'H', 'H'), 'isotopes': (12, 12, 14, 16, 14, 14, 1, 1, 1), 'coords': ((-1.0108, -0.0114, -0.061), (0.478, 0.0191, 0.0139), (1.2974, -0.993, 0.4693), diff --git a/arc/rmgdb_test.py b/arc/rmgdb_test.py index 5a29777f33..3ff298b688 100644 --- a/arc/rmgdb_test.py +++ b/arc/rmgdb_test.py @@ -131,5 +131,6 @@ def test_clean_rmg_database_object(self): self.assertIsNone(self.rmgdb.kinetics) rmgdb.load_rmg_database(rmgdb=self.rmgdb) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))