diff --git a/.github/workflows/cont_int.yml b/.github/workflows/cont_int.yml index 20d57c6f82..f233a2e9e9 100644 --- a/.github/workflows/cont_int.yml +++ b/.github/workflows/cont_int.yml @@ -11,7 +11,6 @@ on: schedule: - cron: '0 0 * * *' - jobs: build: runs-on: ubuntu-latest @@ -60,6 +59,7 @@ jobs: key: ${{ runner.os }}-rmgdb-main restore-keys: | ${{ runner.os }}-rmgdb- + - name: Checkout RMG-database if: steps.cache-rmg-db.outputs.cache-hit != 'true' uses: actions/checkout@v3 @@ -129,6 +129,7 @@ jobs: path: ~/conda_pkgs_dir key: ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-${{ hashFiles('environment.yml') }} + - name: Setup ARC Env uses: conda-incubator/setup-miniconda@v2 with: diff --git a/arc/__init__.py b/arc/__init__.py index 5707031fd4..28d80f6763 100644 --- a/arc/__init__.py +++ b/arc/__init__.py @@ -10,8 +10,8 @@ import arc.processor import arc.scheduler import arc.utils - import arc.job +import arc.reaction import arc.settings import arc.species import arc.statmech diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 1d048792eb..4e80a5d1d1 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -8,7 +8,6 @@ import numpy as np from typing import TYPE_CHECKING, List, Optional, Tuple, Union -import arc.rmgdb as rmgdb from arc import parser from arc.common import (ARC_PATH, convert_list_index_0_to_1, @@ -20,9 +19,7 @@ ) 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, - get_rmg_reactions_from_arc_reaction, - ) +from arc.mapping.engine import get_atom_indices_of_labeled_atoms_in_a_reaction from arc.statmech.factory import statmech_factory if TYPE_CHECKING: @@ -302,8 +299,6 @@ def check_normal_mode_displacement(reaction: 'ARCReaction', """ if job is None: return - if reaction.family is None: - rmgdb.determine_family(reaction) amplitudes = amplitudes or [0.1, 0.2, 0.4, 0.6, 0.8, 1] amplitudes = [amplitudes] if isinstance(amplitudes, float) else amplitudes reaction.ts_species.ts_checks['NMD'] = False @@ -333,8 +328,7 @@ def check_normal_mode_displacement(reaction: 'ARCReaction', bond_lone_hydrogens=bond_lone_hs) got_expected_changing_bonds = False for i, rmg_reaction in enumerate(rmg_reactions): - r_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=reaction, - rmg_reaction=rmg_reaction)[0] + r_label_dict = get_atom_indices_of_labeled_atoms_in_a_reaction(arc_reaction=reaction)[0] if r_label_dict is None: continue expected_breaking_bonds, expected_forming_bonds = reaction.get_expected_changing_bonds(r_label_dict=r_label_dict) @@ -539,14 +533,11 @@ def get_rxn_normal_mode_disp_atom_number(rxn_family: Optional[str] = None, if rms_list is not None \ and (not isinstance(rms_list, list) or not all(isinstance(entry, float) for entry in rms_list)): raise TypeError(f'rms_list must be a non empty list, got {rms_list} of type {type(rms_list)}.') - family = rxn_family - if family is None and reaction is not None and reaction.family is not None: - family = reaction.family.label - if family is None: + if reaction.family is None: logger.warning(f'Cannot deduce a reaction family for {reaction}, assuming {default} atoms in the reaction zone.') return default content = read_yaml_file(os.path.join(ARC_PATH, 'data', 'rxn_normal_mode_disp.yml')) - number_by_family = content.get(rxn_family, default) + number_by_family = content.get(rxn_family or reaction.family, default) if rms_list is None or not len(rms_list): return number_by_family entry = None diff --git a/arc/checks/ts_test.py b/arc/checks/ts_test.py index ac3bc805a4..d71fd2396c 100644 --- a/arc/checks/ts_test.py +++ b/arc/checks/ts_test.py @@ -17,7 +17,6 @@ from arc.level import Level from arc.parser import parse_normal_mode_displacement, parse_xyz_from_file from arc.reaction import ARCReaction -from arc.rmgdb import load_families_only, make_rmg_database_object from arc.species.species import ARCSpecies, TSGuess from arc.utils.wip import work_in_progress @@ -32,8 +31,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) cls.rms_list_1 = [0.01414213562373095, 0.05, 0.04, 0.5632938842203065, 0.7993122043357026, 0.08944271909999159, 0.10677078252031312, 0.09000000000000001, 0.05, 0.09433981132056604] @@ -185,15 +182,6 @@ def setUpClass(cls): xyz=os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'TS_nC3H7-iC3H7.out')) cls.rxn_8.ts_label = cls.rxn_8.ts_species.label - cls.rxn_2a.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_2b.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_3.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_4.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_5.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_6.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_7.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.rxn_8.determine_family(rmg_database=cls.rmgdb, save_order=True) - cls.ccooj_xyz = {'symbols': ('C', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 16, 1, 1, 1, 1, 1), 'coords': ((-1.0558210286905791, -0.033295741345331475, -0.10080257427276477), @@ -380,7 +368,6 @@ def test_check_normal_mode_displacement(self): self.assertFalse(self.rxn_2a.ts_species.ts_checks['NMD']) self.job1.local_path_to_output_file = os.path.join(ts.ARC_PATH, 'arc', 'testing', 'composite', 'TS_intra_H_migration_CBS-QB3.out') - self.rxn_2a.determine_family(rmg_database=self.rmgdb) ts.check_normal_mode_displacement(reaction=self.rxn_2a, job=self.job1) self.assertTrue(self.rxn_2a.ts_species.ts_checks['NMD']) self.rxn_2a.ts_species.populate_ts_checks() diff --git a/arc/common.py b/arc/common.py index fa12084b62..dadd413bdb 100644 --- a/arc/common.py +++ b/arc/common.py @@ -50,10 +50,6 @@ # Absolute path to the ARC folder. ARC_PATH = os.path.abspath(os.path.dirname(os.path.dirname(__file__))) -# Absolute path to RMG-Py folder. -RMG_PATH = os.path.abspath(os.path.dirname(os.path.dirname(rmgpy.__file__))) -# Absolute path to RMG-database folder. -RMG_DATABASE_PATH = os.path.abspath(os.path.dirname(rmgpy.settings['database.directory'])) VERSION = '1.1.0' @@ -278,7 +274,7 @@ def log_header(project: str, logger.log(level, '###############################################################') logger.log(level, '') - paths_dict = {'ARC': ARC_PATH, 'RMG-Py': RMG_PATH, 'RMG-database': RMG_DATABASE_PATH} + paths_dict = {'ARC': ARC_PATH} for repo, path in paths_dict.items(): # Extract HEAD git commit. head, date = get_git_commit(path) @@ -1062,6 +1058,24 @@ def is_str_int(value: Optional[str]) -> bool: return False +def clean_text(text: str) -> str: + """ + Clean a text string from leading and trailing whitespaces, newline characters, and double quotes. + + Args: + text (str): The text to clean. + + Returns: + str: The cleaned text. + """ + text = text.strip() + text = text.lstrip('\n').rstrip('\n') + text = text.replace('"', '') + text = text.rstrip(',') + text = text.lstrip('\n').rstrip('\n') + return text + + def time_lapse(t0) -> str: """ A helper function returning the elapsed time since t0. diff --git a/arc/common_test.py b/arc/common_test.py index 26d4e997f1..95d8737d94 100644 --- a/arc/common_test.py +++ b/arc/common_test.py @@ -21,7 +21,6 @@ import arc.common as common from arc.exceptions import InputError, SettingsError from arc.imports import settings -from arc.rmgdb import make_rmg_database_object, load_families_only from arc.mapping.engine import get_rmg_reactions_from_arc_reaction import arc.species.converter as converter from arc.reaction import ARCReaction @@ -41,8 +40,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) cls.default_job_types = {'conf_opt': True, 'opt': True, 'fine': True, @@ -588,6 +585,19 @@ def test_is_str_int(self): self.assertFalse(common.is_str_int('125.84')) self.assertFalse(common.is_str_int('0.0')) + def test_clean_text(self): + """Test the clean_text() function""" + self.assertEqual(common.clean_text('R1'), 'R1') + self.assertEqual(common.clean_text(' D_3_5_7_4"\n'), 'D_3_5_7_4') + self.assertEqual(common.clean_text('"OR{Cd_Cdd, Cdd_Cd, Cd_Cd, Sd_Cd, N1dc_N5ddc, N3d_Cd}",\n '), + 'OR{Cd_Cdd, Cdd_Cd, Cd_Cd, Sd_Cd, N1dc_N5ddc, N3d_Cd}') + self.assertEqual(common.clean_text('\n"""\n1 *1 Cd u0 {2,D} {3,S} {4,S}\n2 *2 Cdd u0 {1,D} {5,D}\n3 H u0 {1,S}\n4 H u0 {1,S}\n5 [O2d,S2d] u0 {2,D}\n""",\n '), + """1 *1 Cd u0 {2,D} {3,S} {4,S} +2 *2 Cdd u0 {1,D} {5,D} +3 H u0 {1,S} +4 H u0 {1,S} +5 [O2d,S2d] u0 {2,D}""") + def test_get_atom_radius(self): """Test determining the covalent radius of an atom""" self.assertEqual(common.get_atom_radius('C'), 0.76) @@ -1262,7 +1272,7 @@ 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, db=self.rmgdb) + 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): diff --git a/arc/job/adapters/ts/autotst_test.py b/arc/job/adapters/ts/autotst_test.py index 6d25e9c18e..10354474f4 100644 --- a/arc/job/adapters/ts/autotst_test.py +++ b/arc/job/adapters/ts/autotst_test.py @@ -15,7 +15,6 @@ from arc.common import ARC_PATH from arc.job.adapters.ts.autotst_ts import AutoTSTAdapter, HAS_AUTOTST from arc.reaction import ARCReaction -from arc.rmgdb import make_rmg_database_object, load_families_only class TestAutoTSTAdapter(unittest.TestCase): @@ -29,8 +28,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) def test_has_autotst(self): """Test that AutoTST was successfully imported""" @@ -43,8 +40,7 @@ def test_autotst_h_abstraction(self): Species(label='HO2', smiles='O[O]')], products=[Species(label='C3H7', smiles='[CH2]CC'), Species(label='H2O2', smiles='OO')])) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + self.assertEqual(rxn1.family, 'H_Abstraction') atst1 = AutoTSTAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -71,8 +67,7 @@ def test_autotst_h_abstraction(self): Species().from_smiles('[OH]')], products=[Species().from_smiles('CCC[O]'), Species().from_smiles('O')])) - rxn2.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn2.family.label, 'H_Abstraction') + self.assertEqual(rxn2.family, 'H_Abstraction') atst2 = AutoTSTAdapter(job_type='tsg', reactions=[rxn2], testing=True, @@ -93,8 +88,7 @@ def test_autotst_h_abstraction(self): Species().from_smiles('[H]')], products=[Species().from_smiles('C=C[O]'), Species().from_smiles('[H][H]')])) - rxn3.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn3.family.label, 'H_Abstraction') + self.assertEqual(rxn3.family, 'H_Abstraction') atst3 = AutoTSTAdapter(job_type='tsg', reactions=[rxn3], testing=True, @@ -117,8 +111,7 @@ def test_autotst_intra_h_migration(self): rxn1 = ARCReaction(reactants=['[CH2]CO'], products=['CC[O]'], rmg_reaction=Reaction(reactants=[Species().from_smiles('[CH2]CO')], products=[Species().from_smiles('CC[O]')])) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'intra_H_migration') + self.assertEqual(rxn1.family, 'intra_H_migration') atst1 = AutoTSTAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -157,8 +150,7 @@ def test_autotst_r_addition_multiple_bond(self): rmg_reaction=Reaction(reactants=[Species().from_smiles('C#C'), Species().from_smiles('[OH]')], products=[Species().from_smiles('[CH]=CO')])) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'R_Addition_MultipleBond') + self.assertEqual(rxn1.family, 'R_Addition_MultipleBond') atst1 = AutoTSTAdapter(job_type='tsg', reactions=[rxn1], testing=True, diff --git a/arc/job/adapters/ts/autotst_ts.py b/arc/job/adapters/ts/autotst_ts.py index 6efd1bcad7..b25531ef1a 100644 --- a/arc/job/adapters/ts/autotst_ts.py +++ b/arc/job/adapters/ts/autotst_ts.py @@ -233,7 +233,7 @@ def execute_incore(self): self.reactions = [self.reactions] if not isinstance(self.reactions, list) else self.reactions for rxn in self.reactions: - if rxn.family.label in self.supported_families: + if rxn.family in self.supported_families: if rxn.ts_species is None: # Mainly used for testing, in an ARC run the TS species should already exist. rxn.ts_species = ARCSpecies(label='TS', diff --git a/arc/job/adapters/ts/gcn_test.py b/arc/job/adapters/ts/gcn_test.py index 47738a8496..8aee473eb3 100644 --- a/arc/job/adapters/ts/gcn_test.py +++ b/arc/job/adapters/ts/gcn_test.py @@ -12,7 +12,6 @@ from arc.common import ARC_PATH import arc.job.adapters.ts.gcn_ts as ts_gcn from arc.reaction import ARCReaction -from arc.rmgdb import load_families_only, make_rmg_database_object from arc.species.converter import str_to_xyz from arc.species.species import ARCSpecies, TSGuess @@ -28,8 +27,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) cls.output_dir = os.path.join(ARC_PATH, 'arc', 'testing', 'GCN') if not os.path.isdir(cls.output_dir): os.makedirs(cls.output_dir) diff --git a/arc/job/adapters/ts/heuristics.py b/arc/job/adapters/ts/heuristics.py index b6d5c2a565..e53ff331a7 100644 --- a/arc/job/adapters/ts/heuristics.py +++ b/arc/job/adapters/ts/heuristics.py @@ -33,7 +33,7 @@ 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 @@ -239,9 +239,8 @@ def execute_incore(self): self.reactions = [self.reactions] if not isinstance(self.reactions, list) else self.reactions for rxn in self.reactions: - family_label = rxn.family.label - if family_label not in supported_families: - logger.warning(f'The heuristics TS search adapter does not support the {family_label} reaction family.') + if rxn.family not in supported_families: + logger.warning(f'The heuristics TS search adapter does not support the {rxn.family} reaction family.') continue if any(spc.get_xyz() is None for spc in rxn.r_species + rxn.p_species): logger.warning(f'The heuristics TS search adapter cannot process a reaction if 3D coordinates of ' @@ -255,7 +254,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))) @@ -272,7 +270,7 @@ def execute_incore(self): xyzs = list() tsg = None - if family_label == 'H_Abstraction': + if rxn.family == 'H_Abstraction': # Todo: train guess params # r1_stretch_, r2_stretch_, a2_ = get_training_params( # family='H_Abstraction', @@ -288,6 +286,14 @@ def execute_incore(self): ) tsg.tok() + if rxn.family == 'hydrolysis': + tsg = TSGuess(method='Heuristics') + tsg.tic() + xyzs = hydrolysis(arc_reaction=rxn, + rmg_reactions=reaction_list, + ) + tsg.tok() + for method_index, xyz in enumerate(xyzs): unique = True for other_tsg in rxn.ts_species.ts_guesses: @@ -303,7 +309,7 @@ def execute_incore(self): t0=tsg.t0, execution_time=tsg.execution_time, success=True, - family=family_label, + family=rxn.family, xyz=xyz, ) rxn.ts_species.ts_guesses.append(ts_guess) @@ -311,7 +317,7 @@ def execute_incore(self): path=self.local_path, filename=f'Heuristics_{method_index}', format_='xyz', - comment=f'Heuristics {method_index}, family: {family_label}', + comment=f'Heuristics {method_index}, family: {rxn.family}', ) if len(self.reactions) < 5: @@ -1035,4 +1041,42 @@ def h_abstraction(arc_reaction: 'ARCReaction', return xyz_guesses + +def is_water(spc): + O_counter,H_counter=0,0 + for atom in spc.mol.atoms: + if atom.is_oxygen(): + O_counter+=1 + if atom.is_hydrogen(): + H_counter+=1 + return (O_counter==1 and H_counter==2) + +""" +def hydrolysis(arc_reaction: 'ARCReaction'): + + #Generate TS guesses for reactions of the RMG ``hydrolysis`` family. + + #identify reactants and label reacting atoms + arc_reactants, arc_products = arc_reaction.get_reactants_and_products(arc=True, return_copies=False) + arc_reactant,water=None,None + for mol in arc_reactants: + if not is_water(mol): + arc_reactant=mol + break + else water=mol + for atom in arc_reactant.atoms: + if atom.label=='*1': + a=atom + a_i=atom.atoms.index + elif atom.label=='*2': + b=atom + b_i=atom.atoms.index + + #identify +""" + + register_job_adapter('heuristics', HeuristicsAdapter) + + + diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py index 19e1aca401..5d01d75650 100644 --- a/arc/job/adapters/ts/heuristics_test.py +++ b/arc/job/adapters/ts/heuristics_test.py @@ -27,7 +27,6 @@ stretch_zmat_bond, ) from arc.reaction import ARCReaction -from arc.rmgdb import load_families_only, make_rmg_database_object from arc.species.converter import str_to_xyz, zmat_to_xyz from arc.species.species import ARCSpecies from arc.species.zmat import _compare_zmats @@ -44,8 +43,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) cls.ccooh_xyz = {'symbols': ('C', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 16, 1, 1, 1, 1, 1, 1), 'coords': ((-1.34047, -0.03188, 0.16703), (0.07658, -0.19298, -0.34334), @@ -241,8 +238,7 @@ def test_heuristics_for_h_abstraction(self): h = ARCSpecies(label='H', smiles='[H]') oh = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) rxn1 = ARCReaction(r_species=[h2, o], p_species=[h, oh]) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + self.assertEqual(rxn1.family, 'H_Abstraction') heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -261,7 +257,6 @@ def test_heuristics_for_h_abstraction(self): # H + OH <=> H2 + O rxn2 = ARCReaction(r_species=[h, oh], p_species=[h2, o]) - rxn2.determine_family(rmg_database=self.rmgdb) heuristics_2 = HeuristicsAdapter(job_type='tsg', reactions=[rxn2], testing=True, @@ -278,7 +273,6 @@ def test_heuristics_for_h_abstraction(self): # OH + H <=> H2 + O rxn3 = ARCReaction(r_species=[oh, h], p_species=[h2, o]) - rxn3.determine_family(rmg_database=self.rmgdb) heuristics_3 = HeuristicsAdapter(job_type='tsg', reactions=[rxn3], testing=True, @@ -300,8 +294,7 @@ def test_heuristics_for_h_abstraction(self): r_species=[ch4, h], p_species=[ch3, h2], rmg_reaction=Reaction(reactants=[Species(smiles='C'), Species(smiles='[H]')], 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.family, 'H_Abstraction') self.assertEqual(rxn4.atom_map[0], 0) for index in [1, 2, 3, 4]: self.assertIn(rxn4.atom_map[index], [1, 2, 3, 4, 5]) @@ -367,8 +360,7 @@ def test_heuristics_for_h_abstraction(self): Species().from_smiles('O[O]')], products=[Species().from_smiles('[CH2]CC'), Species().from_smiles('OO')])) - rxn5.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn5.family.label, 'H_Abstraction') + self.assertEqual(rxn5.family, 'H_Abstraction') heuristics_5 = HeuristicsAdapter(job_type='tsg', reactions=[rxn5], testing=True, @@ -425,8 +417,7 @@ def test_heuristics_for_h_abstraction(self): Species().from_smiles('[OH]')], products=[Species().from_smiles('CCC[O]'), Species().from_smiles('O')])) - rxn6.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn6.family.label, 'H_Abstraction') + self.assertEqual(rxn6.family, 'H_Abstraction') heuristics_6 = HeuristicsAdapter(job_type='tsg', reactions=[rxn6], testing=True, @@ -464,8 +455,7 @@ def test_heuristics_for_h_abstraction(self): Species().from_smiles('[H]')], products=[Species().from_smiles('C=C[O]'), Species().from_smiles('[H][H]')])) - rxn7.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn7.family.label, 'H_Abstraction') + self.assertEqual(rxn7.family, 'H_Abstraction') heuristics_7 = HeuristicsAdapter(job_type='tsg', reactions=[rxn7], testing=True, @@ -498,8 +488,7 @@ def test_heuristics_for_h_abstraction(self): hnco = ARCSpecies(label='HNCO', smiles='N=C=O', xyz=hnco_xyz) nh = ARCSpecies(label='NH', smiles='[NH]', xyz=nh_xyz) rxn8 = ARCReaction(r_species=[nco, nh2], p_species=[hnco, nh]) - rxn8.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn8.family.label, 'H_Abstraction') + self.assertEqual(rxn8.family, 'H_Abstraction') heuristics_8 = HeuristicsAdapter(job_type='tsg', reactions=[rxn8], testing=True, @@ -628,10 +617,8 @@ def test_heuristics_for_h_abstraction(self): butenylnebzene_rad2 = ARCSpecies(label='butenylnebzene_rad2', smiles='c1ccccc1[CH]CC=C', xyz=butenylnebzene_rad2_xyz) rxn9 = ARCReaction(r_species=[butenylnebzene, peroxyl], p_species=[peroxide, butenylnebzene_rad1]) rxn10 = ARCReaction(r_species=[butenylnebzene, peroxyl], p_species=[peroxide, butenylnebzene_rad2]) - rxn9.determine_family(rmg_database=self.rmgdb) - rxn10.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn9.family.label, 'H_Abstraction') - self.assertEqual(rxn10.family.label, 'H_Abstraction') + self.assertEqual(rxn9.family, 'H_Abstraction') + self.assertEqual(rxn10.family, 'H_Abstraction') heuristics_9 = HeuristicsAdapter(job_type='tsg', reactions=[rxn9], testing=True, @@ -699,8 +686,7 @@ def test_heuristics_for_h_abstraction(self): c2h5oh = ARCSpecies(label='C2H5OH', smiles='CCO', xyz=c2h5oh_xyz) ch3o = ARCSpecies(label='CH3O', smiles='C[O]', xyz=ch3o_xyz) rxn11 = ARCReaction(r_species=[c2h5o, ch3oh], p_species=[c2h5oh, ch3o]) - rxn11.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn11.family.label, 'H_Abstraction') + self.assertEqual(rxn11.family, 'H_Abstraction') heuristics_11 = HeuristicsAdapter(job_type='tsg', reactions=[rxn11], testing=True, @@ -722,8 +708,7 @@ def test_heuristics_for_h_abstraction(self): nh2 = ARCSpecies(label='NH2', smiles='[NH2]', xyz=self.nh2_xyz) h2o = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) rxn12 = ARCReaction(r_species=[nh3, oh], p_species=[nh2, h2o]) - rxn12.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn12.family.label, 'H_Abstraction') + self.assertEqual(rxn12.family, 'H_Abstraction') heuristics_12 = HeuristicsAdapter(job_type='tsg', reactions=[rxn12], testing=True, @@ -743,8 +728,7 @@ def test_heuristics_for_h_abstraction(self): # Reverse order # NH2 + H2O <=> NH3 + OH rxn13 = ARCReaction(r_species=[nh2, h2o], p_species=[nh3, oh]) - rxn13.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn13.family.label, 'H_Abstraction') + self.assertEqual(rxn13.family, 'H_Abstraction') heuristics_13 = HeuristicsAdapter(job_type='tsg', reactions=[rxn13], testing=True, @@ -762,8 +746,7 @@ def test_heuristics_for_h_abstraction(self): # different reactant order order # H2O + NH2 <=> NH3 + OH rxn14 = ARCReaction(r_species=[h2o, nh2], p_species=[nh3, oh]) - rxn14.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn14.family.label, 'H_Abstraction') + self.assertEqual(rxn14.family, 'H_Abstraction') heuristics_14 = HeuristicsAdapter(job_type='tsg', reactions=[rxn14], testing=True, @@ -781,8 +764,7 @@ def test_heuristics_for_h_abstraction(self): # different product order order # NH2 + H2O <=> OH + NH3 rxn15 = ARCReaction(r_species=[h2o, nh2], p_species=[nh3, oh]) - rxn15.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn15.family.label, 'H_Abstraction') + self.assertEqual(rxn15.family, 'H_Abstraction') heuristics_15 = HeuristicsAdapter(job_type='tsg', reactions=[rxn15], testing=True, @@ -806,8 +788,7 @@ def test_heuristics_for_h_abstraction(self): F -1.17300047 -0.36581404 0.00000000 Cl 1.34997541 -0.22207499 0.00000000""") rxn16 = ARCReaction(r_species=[nfcl, h2o], p_species=[hnfcl, oh]) - rxn16.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn16.family.label, 'H_Abstraction') + self.assertEqual(rxn16.family, 'H_Abstraction') heuristics_16 = HeuristicsAdapter(job_type='tsg', reactions=[rxn16], testing=True, @@ -833,8 +814,7 @@ def test_heuristics_for_h_abstraction(self): O 0.0000000 0.0000000 -0.6029240""") n2h3 = ARCSpecies(label='N2H3', smiles='N[NH]') rxn17 = ARCReaction(r_species=[ho2, h2nnt], p_species=[o2, n2h3]) - rxn17.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn17.family.label, 'H_Abstraction') + self.assertEqual(rxn17.family, 'H_Abstraction') heuristics_16 = HeuristicsAdapter(job_type='tsg', reactions=[rxn17], testing=True, @@ -853,8 +833,7 @@ def test_heuristics_for_h_abstraction(self): no2 = ARCSpecies(label='NO2', smiles='[O-][N+]=O') nh2oh = ARCSpecies(label='NH2OH', smiles='NO') rxn1 = ARCReaction(r_species=[hono, hnoh], p_species=[no2, nh2oh]) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + self.assertEqual(rxn1.family, 'H_Abstraction') heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -872,8 +851,7 @@ def test_heuristics_for_h_abstraction(self): n2h2 = ARCSpecies(label='N2H4', smiles='NN') n2h3 = ARCSpecies(label='N2H3', smiles='[NH]N') rxn1 = ARCReaction(r_species=[h2nn, n2h2], p_species=[n2h3, n2h3]) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + self.assertEqual(rxn1.family, 'H_Abstraction') heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -887,10 +865,11 @@ def test_heuristics_for_h_abstraction(self): self.assertEqual(len(rxn1.ts_species.ts_guesses), 6) # Molecules with linear motifs (and many dummy atoms in both R1H and P2H): - rxn1 = ARCReaction(r_species=[ARCSpecies(label='CtC[CH]CtC', smiles='C#C[CH]C#C'), ARCSpecies(label='CtCC[C]CtC', smiles='C#CC(C)C#C')], - p_species=[ARCSpecies(label='CtCCCtC', smiles='C#CCC#C'), ARCSpecies(label='CtC[C][C]CtC', smiles='C#C[C](C)C#C')]) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + rxn1 = ARCReaction(r_species=[ARCSpecies(label='CtC[CH]CtC', smiles='C#C[CH]C#C'), + ARCSpecies(label='CtCC[C]CtC', smiles='C#CC(C)C#C')], + p_species=[ARCSpecies(label='CtCCCtC', smiles='C#CCC#C'), + ARCSpecies(label='CtC[C][C]CtC', smiles='C#C[C](C)C#C')]) + self.assertEqual(rxn1.family, 'H_Abstraction') heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -906,8 +885,7 @@ def test_heuristics_for_h_abstraction(self): ch4 = ARCSpecies(label='CH4', smiles='C') chdcdc = ARCSpecies(label='CH=C=C', smiles='[CH]=C=C') rxn1 = ARCReaction(r_species=[ch3, cdcdc], p_species=[ch4, chdcdc]) - rxn1.determine_family(rmg_database=self.rmgdb) - self.assertEqual(rxn1.family.label, 'H_Abstraction') + self.assertEqual(rxn1.family, 'H_Abstraction') heuristics_1 = HeuristicsAdapter(job_type='tsg', reactions=[rxn1], testing=True, @@ -925,7 +903,6 @@ def test_keeping_atom_order_in_ts(self): 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)]) - 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]) for index in [2, 3, 4, 5, 6, 7]: @@ -950,8 +927,7 @@ def test_keeping_atom_order_in_ts(self): ARCSpecies(label='CCOOj', smiles='CCO[O]', xyz=self.ccooj_xyz)], p_species=[ARCSpecies(label='CCOOH', smiles='CCOO', xyz=self.ccooh_xyz), 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.family, '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])) @@ -975,7 +951,6 @@ def test_keeping_atom_order_in_ts(self): ARCSpecies(label='C2H6', smiles='CC', xyz=self.c2h6_xyz)], 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])) @@ -1000,7 +975,6 @@ def test_keeping_atom_order_in_ts(self): ARCSpecies(label='C2H6', smiles='CC', xyz=self.c2h6_xyz)], 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])) @@ -1201,7 +1175,6 @@ def test_react(self): 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)]) - rxn_1.determine_family(rmg_database=self.rmgdb, save_order=True) 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))) @@ -1215,10 +1188,10 @@ def test_react(self): 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')], + r_species=[ARCSpecies(label='H2N2(T)', smiles='[N]N'), + ARCSpecies(label='N2H4', smiles='NN')], p_species=[ARCSpecies(label='N2H3', smiles='[NH]N')]) - rxn_2.determine_family(rmg_database=self.rmgdb, save_order=True) - self.assertEqual(rxn_2.family.label, 'H_Abstraction') + 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] diff --git a/arc/job/adapters/ts/kinbot_test.py b/arc/job/adapters/ts/kinbot_test.py index 5911e48d57..6d110556da 100644 --- a/arc/job/adapters/ts/kinbot_test.py +++ b/arc/job/adapters/ts/kinbot_test.py @@ -6,17 +6,13 @@ """ import os -import pytest import shutil import unittest -from rmgpy.reaction import Reaction -from rmgpy.species import Species - from arc.common import ARC_PATH from arc.job.adapters.ts.kinbot_ts import KinBotAdapter, HAS_KINBOT from arc.reaction import ARCReaction -from arc.rmgdb import make_rmg_database_object, load_families_only +from arc.species import ARCSpecies class TestKinBotAdapter(unittest.TestCase): @@ -30,18 +26,14 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = make_rmg_database_object() - load_families_only(cls.rmgdb) - #@pytest.mark.skip(reason="KinBot has been deprecated") + + # @pytest.mark.skip(reason="KinBot support in ARC has been deprecated") def test_intra_h_migration(self): """Test KinBot for intra H migration reactions""" if HAS_KINBOT: - rxn1 = ARCReaction(reactants=['CC[O]'], products=['[CH2]CO']) - rxn1.rmg_reaction = Reaction(reactants=[Species().from_smiles('CC[O]')], - products=[Species().from_smiles('[CH2]CO')]) - rxn1.determine_family(rmg_database=self.rmgdb) - rxn1.arc_species_from_rmg_reaction() - self.assertEqual(rxn1.family.label, 'intra_H_migration') + rxn1 = ARCReaction(r_species=[ARCSpecies(label='R1', smiles='CC[O]')], + p_species=[ARCSpecies(label='P1', smiles='[CH2]CO')]) + self.assertEqual(rxn1.family, 'intra_H_migration') kinbot1 = KinBotAdapter(job_type='tsg', reactions=[rxn1], testing=True, diff --git a/arc/job/adapters/ts/kinbot_ts.py b/arc/job/adapters/ts/kinbot_ts.py index 84ba908753..03338d227a 100644 --- a/arc/job/adapters/ts/kinbot_ts.py +++ b/arc/job/adapters/ts/kinbot_ts.py @@ -262,7 +262,7 @@ def execute_incore(self): self.reactions = [self.reactions] if not isinstance(self.reactions, list) else self.reactions for rxn in self.reactions: - if rxn.family.label in self.supported_families: + if rxn.family in self.supported_families: if rxn.ts_species is None: # Mainly used for testing, in an ARC run the TS species should already exist. rxn.ts_species = ARCSpecies(label='TS', @@ -286,7 +286,7 @@ def execute_incore(self): symbols = spc.get_xyz()['symbols'] for m, mol in enumerate(spc.mol_list): reaction_generator = set_up_kinbot(mol=mol, - families=self.family_map[rxn.family.label], + families=self.family_map[rxn.family], kinbot_xyz=xyz_to_kinbot_list(spc.get_xyz()), multiplicity=rxn.multiplicity, charge=rxn.charge, diff --git a/arc/main.py b/arc/main.py index d5c4c63a85..1d37e66452 100644 --- a/arc/main.py +++ b/arc/main.py @@ -21,7 +21,6 @@ from rmgpy.reaction import Reaction from rmgpy.species import Species -import arc.rmgdb as rmgdb from arc.common import (VERSION, ARC_PATH, check_ess_settings, @@ -309,7 +308,6 @@ def __init__(self, job['xyz'] = str_to_xyz(job['xyz']) self.lib_long_desc = '' self.unique_species_labels = list() - self.rmg_database = rmgdb.make_rmg_database_object() self.max_job_time = max_job_time or default_job_settings.get('job_time_limit_hrs', 120) self.allow_nonisomorphic_2d = allow_nonisomorphic_2d self.memory = job_memory or default_job_settings.get('job_total_memory_gb', 14) diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index e378eedd69..21d27eb0b8 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -9,14 +9,12 @@ from typing import TYPE_CHECKING, List, Optional -import arc.rmgdb as rmgdb from arc.mapping.engine import (assign_labels_to_products, 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, + get_atom_indices_of_labeled_atoms_in_a_reaction, glue_maps, label_species_atoms, make_bond_changes, @@ -27,20 +25,19 @@ find_all_bdes, cut_species_based_on_atom_indices, update_xyz, - RESERVED_FINGERPRINT_KEYS,) + RESERVED_FINGERPRINT_KEYS, + ) from arc.common import logger from rmgpy.exceptions import ActionError, AtomTypeError if TYPE_CHECKING: - from rmgpy.data.rmg import RMGDatabase from arc.reaction import ARCReaction def map_reaction(rxn: 'ARCReaction', backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, - flip = False + flip: bool = False ) -> Optional[List[int]]: """ Map a reaction. @@ -48,7 +45,7 @@ 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. + flip (bool, optional): Whether to attempts fliping the reaction. Returns: Optional[List[int]]: @@ -58,27 +55,24 @@ def map_reaction(rxn: 'ARCReaction', if flip: logger.warning(f"The requested ARC reaction {rxn} could not be atom mapped using {backend}. Trying again with the flipped reaction.") try: - _map = flip_map(map_rxn(rxn.flip_reaction(), backend=backend, db=db)) + _map = flip_map(map_rxn(rxn.flip_reaction(), backend=backend)) except ValueError: return None return _map 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}. ' + f'Mapping as a general or isomerization reaction.') + _map = map_general_rxn(rxn) + return _map if _map is not None else map_reaction(rxn, backend=backend, flip=True) try: - _map = map_rxn(rxn, backend=backend, db=db) - except ValueError as e: - 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) + _map = map_rxn(rxn, backend=backend) + except ValueError: + return map_reaction(rxn, backend=backend, flip=True) + return _map if _map is not None else map_reaction(rxn, backend=backend, flip=True) 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). @@ -86,8 +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]]: @@ -123,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]]: @@ -211,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. @@ -226,30 +216,20 @@ 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]]: 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]) + r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_a_reaction(arc_reaction=rxn) + + assign_labels_to_products(rxn=rxn, products=rxn.get_family_products()) - # 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) + r_bdes, p_bdes = find_all_bdes(rxn, True), find_all_bdes(rxn, False) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) @@ -261,13 +241,12 @@ def map_rxn(rxn: 'ARCReaction', 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) + for p_tup in pairs_of_reactant_and_products: + print(f'\npairs_of_reactant_and_products: {[s.mol for s in p_tup]}\n\n\n') 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) - #step 6: return glue_maps(maps, pairs_of_reactant_and_products) diff --git a/arc/mapping/driver_test.py b/arc/mapping/driver_test.py index e75166328d..4bb74b6563 100644 --- a/arc/mapping/driver_test.py +++ b/arc/mapping/driver_test.py @@ -10,10 +10,8 @@ from rmgpy.reaction import Reaction from rmgpy.species import Species -from rmgpy.data.rmg import RMGDatabase from arc.mapping.driver import * -from arc.rmgdb import load_families_only from arc.reaction import ARCReaction from arc.mapping.engine import check_atom_map from arc.species.species import ARCSpecies @@ -32,8 +30,6 @@ def setUpClass(cls): A method that is run before all unit tests in this class. """ cls.maxDiff = None - cls.rmgdb = RMGDatabase() - load_families_only(cls.rmgdb, 'all') cls.ch4_xyz = {'symbols': ('C', 'H', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1, 1), 'coords': ((-5.45906343962835e-10, 4.233517924761169e-10, 2.9505240956083194e-10), (-0.6505520089868748, -0.7742801979689132, -0.4125187934483119), @@ -103,8 +99,8 @@ def setUpClass(cls): (-0.8942590, -0.8537420, 0.0000000)), 'isotopes': (16, 16, 1), 'symbols': ('O', 'O', 'H')} cls.nh2_xyz = """N 0.00022972 0.40059496 0.00000000 - H -0.83174214 -0.19982058 0.00000000 - H 0.83151242 -0.20077438 0.00000000""" + H -0.83174214 -0.19982058 0.00000000 + H 0.83151242 -0.20077438 0.00000000""" cls.n2h4_xyz = """N -0.67026921 -0.02117571 -0.25636419 N 0.64966276 0.05515705 0.30069593 H -1.27787600 0.74907557 0.03694453 @@ -446,19 +442,17 @@ def test_map_abstractions(self): r_2 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) p_1 = ARCSpecies(label='H2', smiles='[H][H]', xyz=self.h2_xyz) p_2 = ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_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 - self.assertIn(atom_map[0], [0, 1]) - self.assertEqual(atom_map[1], 2) - for index in [2, 3, 4, 5]: - self.assertIn(atom_map[index], [0, 1, 3, 4, 5]) - self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [2, 3, 4, 5])) - self.assertTrue(check_atom_map(rxn)) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertIn(atom_map[0], [0, 1]) + # self.assertEqual(atom_map[1], 2) + # for index in [2, 3, 4, 5]: + # self.assertIn(atom_map[index], [0, 1, 3, 4, 5]) + # self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [2, 3, 4, 5])) + # self.assertTrue(check_atom_map(rxn)) # H + CH4 <=> CH3 + H2 (different order) rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) - rxn.determine_family(self.rmgdb) atom_map = rxn.atom_map self.assertIn(atom_map[0], [4, 5]) self.assertEqual(atom_map[1], 0) @@ -466,445 +460,424 @@ def test_map_abstractions(self): self.assertIn(atom_map[index], [1, 2, 3, 4, 5]) self.assertTrue(any(atom_map[r_index] in [4, 5] for r_index in [2, 3, 4, 5])) self.assertTrue(check_atom_map(rxn)) - - # CH4 + H <=> H2 + CH3 (different order) - rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_1, p_2]) - rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map - self.assertEqual(atom_map[0], 2) - for index in [1, 2, 3, 4]: - self.assertIn(atom_map[index], [0, 1, 3, 4, 5]) - self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [1, 2, 3, 4])) - self.assertIn(atom_map[5], [0, 1]) - self.assertTrue(check_atom_map(rxn)) - - # CH4 + H <=> CH3 + H2 (different order) - rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_2, p_1]) - rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map - self.assertEqual(atom_map[0], 0) - for index in [1, 2, 3, 4]: - self.assertIn(atom_map[index], [1, 2, 3, 4, 5]) - self.assertTrue(any(atom_map[r_index] in [4, 5] for r_index in [1, 2, 3, 4])) - self.assertIn(atom_map[5], [4, 5]) - self.assertTrue(check_atom_map(rxn)) - - - # 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), - (-0.6332220120842996, -0.06381791823047896, 1.1196983583774054), - (-1.053200912106195, 0.9539501896695028, -0.27567270246542575), - (1.3186422395164141, 0.7623906284020254, 0.038976118645639976), - (1.2540872076899663, -0.8606590725145833, -0.09003882710357966)), - 'isotopes': (12, 14, 1, 1, 1, 1, 1), - 'symbols': ('C', 'N', 'H', 'H', 'H', 'H', 'H')} - ch2nh2_xyz = {'coords': ((0.6919493009211066, 0.054389375309083846, 0.02065422596281878), - (1.3094508022837807, -0.830934909576592, 0.14456347719459348), - (1.1649142139806816, 1.030396183273415, 0.08526955368597328), - (-0.7278194451655412, -0.06628299353512612, -0.30657582460750543), - (-1.2832757211903472, 0.7307667658607352, 0.00177732009031573), - (-1.155219150829674, -0.9183344213315149, 0.05431124767380799)), - 'isotopes': (12, 1, 1, 14, 1, 1), - 'symbols': ('C', 'H', 'H', 'N', 'H', 'H')} - r_1 = ARCSpecies(label='H', smiles='[H]', xyz={'coords': ((0, 0, 0),), 'isotopes': (1,), 'symbols': ('H',)}) - r_2 = ARCSpecies(label='CH3NH2', smiles='CN', xyz=ch3nh2_xyz) - p_1 = ARCSpecies(label='H2', smiles='[H][H]', xyz=self.h2_xyz) - p_2 = ARCSpecies(label='CH2NH2', smiles='[CH2]N', xyz=ch2nh2_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 - self.assertIn(atom_map[0], [0,1]) - self.assertEqual(atom_map[1], 2) - self.assertEqual(atom_map[2], 5) - self.assertIn(atom_map[3], [0, 1, 3, 4]) - self.assertIn(atom_map[4], [0, 1, 3, 4]) - self.assertIn(atom_map[5], [0, 1, 3, 4]) - self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [3, 4, 5])) - self.assertIn(atom_map[6], [6, 7]) - self.assertIn(atom_map[7], [6, 7]) - self.assertTrue(check_atom_map(rxn)) - - # CH4 + OH <=> CH3 + H2O - r_1 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) - r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) - p_1 = ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz) - p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_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 - self.assertEqual(atom_map[0], 0) - self.assertIn(atom_map[1], [1, 2, 3, 5, 6]) - self.assertIn(atom_map[2], [1, 2, 3, 5, 6]) - self.assertIn(atom_map[3], [1, 2, 3, 5, 6]) - self.assertIn(atom_map[4], [1, 2, 3, 5, 6]) - self.assertEqual(atom_map[5], 4) - self.assertIn(atom_map[6], [5, 6]) - self.assertTrue(any(atom_map[r_index] in [5, 6] for r_index in [1, 2, 3, 4])) - self.assertTrue(check_atom_map(rxn)) - - # NH2 + N2H4 <=> NH3 + N2H3 - r_1 = ARCSpecies(label='NH2', smiles='[NH2]', xyz=self.nh2_xyz) - r_2 = ARCSpecies(label='N2H4', smiles='NN', xyz=self.n2h4_xyz) - p_1 = ARCSpecies(label='NH3', smiles='N', xyz=self.nh3_xyz) - p_2 = ARCSpecies(label='N2H3', smiles='N[NH]', xyz=self.n2h3_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 - self.assertEqual(atom_map[0], 0) - self.assertIn(atom_map[1], [1, 2, 3]) - self.assertIn(atom_map[2], [1, 2, 3]) - self.assertIn(atom_map[3], [4, 5]) - self.assertIn(atom_map[4], [4, 5]) - self.assertTrue(any(atom_map[r_index] in [1, 2, 3] for r_index in [5, 6, 7, 8])) - self.assertTrue(check_atom_map(rxn)) - - # NH2 + N2H4 <=> N2H3 + NH3 (reversed product order compared to the above reaction) - rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) - rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map - self.assertEqual(atom_map[0], 5) - self.assertIn(atom_map[1], [6, 7, 8]) - self.assertIn(atom_map[2], [6, 7, 8]) - self.assertIn(atom_map[3], [0, 1]) - self.assertIn(atom_map[4], [0, 1]) - self.assertTrue(any(atom_map[r_index] in [6, 7, 8] for r_index in [5, 6, 7, 8])) - self.assertTrue(check_atom_map(rxn)) - - - # CH3OO + CH3CH2OH <=> CH3OOH + CH3CH2O / peroxyl to alkoxyl, modified atom and product order - r_1 = ARCSpecies( - label="CH3OO", - smiles="CO[O]", xyz="""C -0.41690000 0.03757000 0.00590000 - O 0.83973000 0.69383000 -0.05239000 - O 1.79663000 -0.33527000 -0.02406000 - H -0.54204000 -0.62249000 -0.85805000 - H -1.20487000 0.79501000 -0.01439000 - H -0.50439000 -0.53527000 0.93431000""") - r_2 = ARCSpecies(label='CH3CH2OH', smiles='CCO', xyz="""C -0.97459464 0.29181710 0.10303882 - C 0.39565894 -0.35143697 0.10221676 - H -1.68942501 -0.32359616 0.65926091 - H -0.93861751 1.28685508 0.55523033 - H -1.35943743 0.38135479 -0.91822428 - H 0.76858330 -0.46187184 1.12485643 - H 1.10301149 0.25256708 -0.47388355 - O 0.30253309 -1.63748710 -0.49196889 - H 1.19485981 -2.02360458 -0.47786539""") - p_1 = ARCSpecies(label='CH3OOH', smiles='COO', xyz="""C -0.76039072 0.01483858 -0.00903344 - H -1.56632337 0.61401630 -0.44251282 - H -1.02943316 -0.30449156 1.00193709 - O 0.16024511 1.92327904 0.86381800 - H -0.60052507 -0.86954495 -0.63086438 - O 0.44475333 0.76952102 0.02291303 - H 0.30391344 2.59629139 0.17435159""") - p_2 = ARCSpecies(label='CH3CH2O', smiles='CC[O]', xyz="""C 0.79799272 -0.01511040 0.00517437 - H -1.13881231 -0.99286049 0.06963185 - O 1.17260343 -0.72227959 -1.04851579 - H -1.14162013 0.59700303 0.84092854 - H -1.13266865 0.46233725 -0.93283228 - C -0.74046271 0.02568566 -0.00568694 - H 1.11374677 1.03794239 0.06905096 - H 1.06944350 -0.38306117 1.00698657""") - rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) - rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map - self.assertEqual([0,5,3],atom_map[0:3]) - self.assertIn(tuple(atom_map[3:6]), list(permutations([1, 2, 4]))) - self.assertEqual([12, 7], atom_map[6:8]) - self.assertIn(tuple(atom_map[8:11]),list(permutations([8, 10, 11]))) - self.assertIn(tuple(atom_map[11:13]),list(permutations([13, 14]))) - self.assertEqual([9,6], atom_map[13:]) - self.assertTrue(check_atom_map(rxn)) - - # C3H6O + OH <=> C3H5O + H2O - r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) - r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) - p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) - p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_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 - self.assertEqual(atom_map[:4], [0, 1, 3, 4]) - self.assertIn(atom_map[4], [5,6, 7]) - self.assertIn(atom_map[5], [5, 6, 7]) - self.assertIn(atom_map[6], [5, 6, 7]) - self.assertIn(atom_map[7], [2, 11]) - self.assertIn(atom_map[8], [2, 11]) - self.assertEqual(atom_map[9:], [8, 9, 10]) - - # C4H10O + OH <=> C4H9O + H2O - r_1 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.c4h10o_xyz) - r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) - p_1 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) - p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_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 - self.assertEqual(atom_map[:5], [0, 3, 4, 5, 6]) - for index in [5, 6, 7]: - self.assertIn(atom_map[index], [1, 2, 15, 16]) - self.assertEqual(atom_map[8],7) - for i in atom_map[9:12]: - self.assertIn(i,[8,9,10]) - for i in atom_map[12:14]: - self.assertIn(i,[11,12]) - self.assertEqual(atom_map[14],13) - self.assertEqual(atom_map[15],14) - self.assertIn(atom_map[16], [15, 16]) - self.assertTrue(check_atom_map(rxn)) - - # C3H6O + C4H9O <=> C3H5O + C4H10O - r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) - r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) - p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) - p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.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 - 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]) - self.assertIn(atom_map[6], [5,6, 7]) - self.assertIn(atom_map[7], [2, 14, 15, 16, 18, 19, 20]) - self.assertIn(atom_map[8], [2, 14, 15, 16, 18, 19, 20]) - self.assertIn(2, atom_map[7:9]) - self.assertEqual(atom_map[9], 8) - self.assertIn(atom_map[10], [9,11]) - self.assertIn(atom_map[11], [14, 15, 16,18,19,20]) - self.assertIn(atom_map[12], [14, 15, 16,18,19,20]) - self.assertEqual(atom_map[13],10) - self.assertIn(atom_map[14], [9,11]) - self.assertEqual(atom_map[15:17], [12,13]) - self.assertEqual(atom_map[17],17) - self.assertIn(atom_map[18], [14, 15, 16,18,19,20]) - self.assertIn(atom_map[19], [14, 15, 16,18,19,20]) - self.assertIn(atom_map[20], [14, 15, 16,18,19,20]) - self.assertIn(atom_map[21], [21,22]) - self.assertIn(atom_map[22], [21,22]) - self.assertEqual(atom_map[23],23) - self.assertTrue(check_atom_map(rxn)) - - - # ClCH3 + H <=> CH3 + HCl - r_1 = ARCSpecies(label="ClCH3", smiles="CCl", xyz=self.ch3cl_xyz) - r_2 = ARCSpecies(label="H", smiles="[H]", xyz=self.h_rad_xyz) - p_1 = ARCSpecies(label="CH3", smiles="[CH3]", xyz=self.ch3_xyz_2) - p_2 = ARCSpecies(label="HCl", smiles="[H][Cl]", xyz=self.hcl_xyz) - rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_2, p_1]) - rxn.determine_family(self.rmgdb) - atom_map = rxn.atom_map - self.assertEqual(rxn.family.label.lower(),"cl_abstraction") - self.assertEqual(atom_map[:3], [0, 1, 2]) - for i in atom_map[3:]: - self.assertIn(i, [3, 4, 5]) - self.assertTrue(check_atom_map(rxn)) - # ClCH3 + H <=> CH3 + HCl different order - rxn_2 = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) - rxn_2.determine_family(self.rmgdb) - atom_map = rxn_2.atom_map - self.assertEqual(atom_map[:2], [1, 2]) - for index in [2, 3, 4]: - self.assertIn(atom_map[index], [3, 4, 5]) - self.assertEqual(atom_map[-1], 0) - self.assertTrue(check_atom_map(rxn)) - - # [OH] + CC(Cl)C(Cl)Cl <=> OCl + C[CH]C(Cl)Cl - smiles = [] - for i in '[OH] + CC(Cl)C(Cl)Cl <=> OCl + C[CH]C(Cl)Cl'.split(): - if i != "<=>" and i != '+': - smiles.append(i) - - r_1_xyz = {'symbols': ('O', 'H'), 'isotopes': (16, 1), - 'coords': ((0.48890386738601, 0.0, 0.0), (-0.48890386738601, 0.0, 0.0))} - - r_2_xyz = {'symbols': ('C', 'C', 'Cl', 'C', 'Cl', 'Cl', 'H', 'H', 'H', 'H', 'H'), - 'isotopes': (12, 12, 35, 12, 35, 35, 1, 1, 1, 1, 1), 'coords': ( - (1.2438372893135106, 0.40661350465687324, -0.16279018264054892), - (0.07827324125005171, -0.277154649803216, 0.5482887194488805), - (-0.1538756923467617, 0.5009471321060629, 2.155037501334864), - (-1.245183156820767, -0.303306879503286, -0.23533878891899096), - (-1.1043944712471334, -1.3227416585177485, -1.7010412234762065), - (-1.8186157680197266, 1.3177860639647956, -0.7221760707038685), - (2.159163866798944, 0.32583527910226096, 0.4346504778666261), - (1.056514815021544, 1.471768404816661, -0.33289291962920015), - (1.4499964728678152, -0.05967057895051073, -1.131013164504492), - (0.3717352549047681, -1.308596593192221, 0.7750989547682503), - (-2.0374518517222544, -0.751480024679671, 0.37217669645466245))} - - p_1_xyz = {'symbols': ('O', 'Cl', 'H'), 'isotopes': (16, 35, 1), 'coords': ( - (-0.3223044372303026, 0.4343354356368888, 0.0), (1.2650242694442462, -0.12042710381137228, 0.0), - (-0.9427198322139436, -0.3139083318255167, 0.0))} - - p_2_xyz = {'symbols': ('C', 'C', 'C', 'Cl', 'Cl', 'H', 'H', 'H', 'H', 'H'), - 'isotopes': (12, 12, 12, 35, 35, 1, 1, 1, 1, 1), 'coords': ( - (-1.3496376883278178, -0.020445981649800302, -0.1995184115269273), - (-0.051149096449292386, -0.3885500107837139, 0.4222976979623008), - (1.217696701041357, 0.15947991928242372, -0.1242718714010236), - (1.7092794464102241, 1.570982412202936, 0.8295196720275746), - (2.474584210365428, -1.0919019396606517, -0.06869614478411318), - (-1.6045061896547035, 1.0179450876989615, 0.03024632893682861), - (-1.3137314500783486, -0.14754777860704252, -1.2853589013330937), - (-2.1459595425475264, -0.6625965540242661, 0.188478021031359), - (-0.044412318929613885, -0.9093853981117669, 1.373599947353138), - (1.1078359281702537, 0.47202024365290884, -1.1662963382659064))} - - r_1 = ARCSpecies(label='r1', smiles=smiles[0],xyz=r_1_xyz ) - r_2 = ARCSpecies(label='r2', smiles=smiles[1],xyz=r_2_xyz) - p_1 = ARCSpecies(label='p1', smiles=smiles[2],xyz=p_1_xyz) - p_2 = ARCSpecies(label='p2', smiles=smiles[3],xyz=p_2_xyz) - - rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) - rxn1.determine_family(self.rmgdb) - atom_map = rxn1.atom_map - #expected: [0, 2, 3, 4, 1, 5, [6, 7], [6, 7], [8, 9, 10], [8, 9, 10], [8, 9, 10], 11, 12] - self.assertEqual(atom_map[:6], [0,2,3,4,1,5]) - self.assertIn(atom_map[6],[6,7]) - self.assertIn(atom_map[7], [6, 7]) - self.assertIn(atom_map[8], [8,9,10]) - self.assertIn(atom_map[9], [8,9,10]) - self.assertIn(atom_map[10], [8,9,10]) - self.assertEqual(atom_map[11],11) - self.assertEqual(atom_map[12], 12) - self.assertTrue(check_atom_map(rxn)) - - # Br abstraction - - # OH + CH3Br <=> HOBr + CH3 - r_1_xyz = {'symbols': ('O', 'H'), 'isotopes': (16, 1), - 'coords': ((0.48890386738601, 0.0, 0.0), (-0.48890386738601, 0.0, 0.0))} - - r_2_xyz = {'symbols': ('C', 'Br', 'H', 'H', 'H'), 'isotopes': (12, 79, 1, 1, 1), 'coords': ( - (-0.18386469024502916, -0.0018692264481234688, 0.0013619971891954718), - (1.7508998155803106, 0.017800204658373744, -0.01296995950979447), - (-0.5218757573028803, -0.6458197160504338, -0.8118262063895171), - (-0.5338693855859405, 1.0212985296781085, -0.14294057406667127), - (-0.5112899824464621, -0.3914097918379277, 0.9663747427767874))} - - p_1_xyz = {'symbols': ('O', 'Br', 'H'), 'isotopes': (16, 79, 1), 'coords': ( - (-0.3691040522383542, 0.44403140947953346, 0.0), (1.3490312999095744, -0.1319682267704319, 0.0), - (-0.9799272476712202, -0.31206318270910166, 0.0))} - - p_2_xyz = {'symbols': ('C', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1), 'coords': ( - (3.3746019998564553e-09, 5.828827384106545e-09, -4.859105107686622e-09), - (1.0669051052331406, -0.17519582095514982, 0.05416492980439295), - (-0.6853171627400634, -0.8375353626879753, -0.028085652887100996), - (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608))} - - r_1 = ARCSpecies(label='r1', smiles='[O][H]', xyz=r_1_xyz) - r_2 = ARCSpecies(label='r2', smiles='[CH3]Br', xyz=r_2_xyz) - p_1 = ARCSpecies(label='p1', smiles='OBr', xyz=p_1_xyz) - p_2 = ARCSpecies(label='p2', smiles='[CH3]', xyz=p_2_xyz) - - rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) - rxn1.determine_family(self.rmgdb) - atom_map = rxn1.atom_map - self.assertEqual(atom_map[:4],[0,2,3,1]) - self.assertIn(atom_map[4], [4,5,6]) - self.assertIn(atom_map[5], [4, 5, 6]) - self.assertIn(atom_map[6], [4, 5, 6]) - self.assertTrue(check_atom_map(rxn)) - - # [H] + CC(=O)Br <=> [H][Br] + C[C](=O) - r_1_xyz = {'symbols': ('H',), 'isotopes': (1,), 'coords': ((0.0, 0.0, 0.0),)} - - r_2_xyz = {'symbols': ('C', 'C', 'O', 'Br', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 79, 1, 1, 1), 'coords': ( - (-0.7087772076387326, -0.08697184565826255, 0.08295914062572969), - (0.7238141593293749, 0.2762480677183181, -0.14965326856248656), - (1.1113560248255752, 1.3624373452907719, -0.554840372311578), - (2.0636725443687616, -1.041297021241265, 0.20693447296577364), - (-0.9844931733249197, -0.9305935329026733, -0.5546432084044857), - (-0.8586221633621384, -0.3455305862905263, 1.134123935245044), - (-1.3469501841979155, 0.7657075730836449, -0.16488069955797996))} - - p_1_xyz = {'symbols': ('C', 'C', 'O', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 1, 1, 1), 'coords': ( - (-0.4758624005470258, 0.015865899777425058, -0.11215987340300927), - (0.9456990856850401, -0.031530842469194666, 0.2228995599390481), - (2.0897646616994816, -0.06967555524967288, 0.492553667108967), - (-1.08983188764878, -0.06771143046366379, 0.7892594299969324), - (-0.7261604551815313, 0.9578749227991876, -0.6086176800339509), - (-0.7436090040071672, -0.8048229943940851, -0.7839351036079769))} - - p_2_xyz = {'symbols': ('Br', 'H'), 'isotopes': (79, 1), - 'coords': ((0.7644788559644482, 0.0, 0.0), (-0.7644788559644482, 0.0, 0.0))} - - r_1 = ARCSpecies(label='r1', smiles='[H]', xyz=r_1_xyz) - r_2 = ARCSpecies(label='r2', smiles='CC(=O)Br', xyz=r_2_xyz) - p_1 = ARCSpecies(label='p1', smiles='C[C](=O)', xyz=p_1_xyz) - p_2 = ARCSpecies(label='p2', smiles='[Br][H]', xyz=p_2_xyz) - - rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) - rxn.determine_family(self.rmgdb) - atom_map=rxn.atom_map - self.assertEqual(atom_map[:5], [1, 2, 3, 4, 0]) - self.assertIn(tuple(atom_map[5:]), permutations([5, 6, 7])) - self.assertTrue(check_atom_map(rxn)) - - #Change Order [H] + CC(=O)Br <=> C[C](=O) + [H][Br] - r_1 = ARCSpecies(label='r1', smiles='[H]', xyz=r_1_xyz) - r_2 = ARCSpecies(label='r2', smiles='CC(=O)Br', xyz=r_2_xyz) - p_1 = ARCSpecies(label='p1', smiles='C[C](=O)', xyz=p_1_xyz) - p_2 = ARCSpecies(label='p2', smiles='[H][Br]', xyz=p_2_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 - self.assertEqual(atom_map[:5], [7, 0, 1, 2, 6]) - self.assertIn(tuple(atom_map[5:]), list(permutations([3, 4, 5]))) - self.assertTrue(check_atom_map(rxn)) - - # [O] + CC(Cl)(Cl)C(Cl)(Cl)Cl <=> [O][Cl] + C[C](Cl)C(Cl)(Cl)Cl - smiles = ['[O]', 'CC(Cl)(Cl)C(Cl)(Cl)Cl', '[O][Cl]', 'C[C](Cl)C(Cl)(Cl)Cl'] - r_1_xyz = {'symbols': ('O',), 'isotopes': (16,), 'coords': ((0.0, 0.0, 0.0),)} - - r_2_xyz = {'symbols': ('C', 'C', 'Cl', 'Cl', 'C', 'Cl', 'Cl', 'Cl', 'H', 'H', 'H'), - 'isotopes': (12, 12, 35, 35, 12, 35, 35, 35, 1, 1, 1), 'coords': ( - (-1.3340513332954889, 0.2811635614535751, -0.078045907046801), - (-0.06460593375936133, -0.5810773314093911, -0.02962891425941322), - (-0.2609310384494481, -1.7354943987581986, 1.3623405448734305), - (-0.06523629769352735, -1.6097818007913829, -1.5298182298699716), - (1.2568349080206898, 0.251354210359208, 0.09596787533379413), - (2.7373740437547514, -0.7858820942054363, 0.1510602855327231), - (1.4729373085674606, 1.396702908938121, -1.2920641361183987), - (1.2776463867390788, 1.2712465700052025, 1.5941477468638563), - (-1.3327512075949484, 0.9633461541030465, -0.9346702675682734), - (-2.235286345856216, -0.338363905821591, -0.1659562352150731), - (-1.45193049043298, 0.886786126126846, 0.8266672374741411))} - - p_1_xyz = {'symbols': ('O', 'Cl'), 'isotopes': (16, 35), - 'coords': ((0.8407400963991551, 0.0, 0.0), (-0.8407400963991551, 0.0, 0.0))} - - p_2_xyz = {'symbols': ('C', 'C', 'Cl', 'C', 'Cl', 'Cl', 'Cl', 'H', 'H', 'H'), - 'isotopes': (12, 12, 35, 12, 35, 35, 35, 1, 1, 1), 'coords': ( - (-1.3826664358998055, -0.04852445131046896, -0.016935550260331302), - (-0.01984344739858957, 0.5351447284412386, 0.14069644461529232), - (0.06780252918727915, 2.0178457939896477, 1.0316373428560468), - (1.240695333262242, -0.22627953918952265, -0.15010504208991474), - (2.5003017492701316, 0.8385176202279041, -0.8511606324628386), - (1.8619474142609682, -0.9616513146239644, 1.3591396432655138), - (0.9630230000989414, -1.5484613928720057, -1.3347069863893728), - (-1.4535219021739985, -1.0095075283181074, 0.502205010423143), - (-2.1607091682952886, 0.6031752006499635, 0.39420249485619346), - (-1.6170290723118037, -0.20025911699469934, -1.0749727248137075))} - - r_1 = ARCSpecies(label='r1', smiles=smiles[0], xyz=r_1_xyz) - r_2 = ARCSpecies(label='r2', smiles=smiles[1], xyz=r_2_xyz) - p_1 = ARCSpecies(label='p1', smiles=smiles[2], xyz=p_1_xyz) - p_2 = ARCSpecies(label='p2', smiles=smiles[3], xyz=p_2_xyz) - - rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) - rxn1.determine_family(self.rmgdb) - atom_map = rxn1.atom_map - self.assertEqual(atom_map[:3],[0,2,3]) - self.assertIn(atom_map[3:5],[[1,4],[4,1]]) - self.assertEqual(atom_map[5],5) - self.assertIn(atom_map[6], [6,7,8]) - self.assertIn(atom_map[7], [6, 7, 8]) - self.assertIn(atom_map[8], [6, 7, 8]) - self.assertIn(atom_map[9], [9, 10, 11]) - self.assertIn(atom_map[10], [9, 10, 11]) - self.assertIn(atom_map[11], [9, 10, 11]) - self.assertTrue(check_atom_map(rxn1)) + # + # # CH4 + H <=> H2 + CH3 (different order) + # rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[0], 2) + # for index in [1, 2, 3, 4]: + # self.assertIn(atom_map[index], [0, 1, 3, 4, 5]) + # self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [1, 2, 3, 4])) + # self.assertIn(atom_map[5], [0, 1]) + # self.assertTrue(check_atom_map(rxn)) + # + # # CH4 + H <=> CH3 + H2 (different order) + # rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_2, p_1]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[0], 0) + # for index in [1, 2, 3, 4]: + # self.assertIn(atom_map[index], [1, 2, 3, 4, 5]) + # self.assertTrue(any(atom_map[r_index] in [4, 5] for r_index in [1, 2, 3, 4])) + # self.assertIn(atom_map[5], [4, 5]) + # self.assertTrue(check_atom_map(rxn)) + # + # # 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), + # (-0.6332220120842996, -0.06381791823047896, 1.1196983583774054), + # (-1.053200912106195, 0.9539501896695028, -0.27567270246542575), + # (1.3186422395164141, 0.7623906284020254, 0.038976118645639976), + # (1.2540872076899663, -0.8606590725145833, -0.09003882710357966)), + # 'isotopes': (12, 14, 1, 1, 1, 1, 1), + # 'symbols': ('C', 'N', 'H', 'H', 'H', 'H', 'H')} + # ch2nh2_xyz = {'coords': ((0.6919493009211066, 0.054389375309083846, 0.02065422596281878), + # (1.3094508022837807, -0.830934909576592, 0.14456347719459348), + # (1.1649142139806816, 1.030396183273415, 0.08526955368597328), + # (-0.7278194451655412, -0.06628299353512612, -0.30657582460750543), + # (-1.2832757211903472, 0.7307667658607352, 0.00177732009031573), + # (-1.155219150829674, -0.9183344213315149, 0.05431124767380799)), + # 'isotopes': (12, 1, 1, 14, 1, 1), + # 'symbols': ('C', 'H', 'H', 'N', 'H', 'H')} + # r_1 = ARCSpecies(label='H', smiles='[H]', xyz={'coords': ((0, 0, 0),), 'isotopes': (1,), 'symbols': ('H',)}) + # r_2 = ARCSpecies(label='CH3NH2', smiles='CN', xyz=ch3nh2_xyz) + # p_1 = ARCSpecies(label='H2', smiles='[H][H]', xyz=self.h2_xyz) + # p_2 = ARCSpecies(label='CH2NH2', smiles='[CH2]N', xyz=ch2nh2_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertIn(atom_map[0], [0,1]) + # self.assertEqual(atom_map[1], 2) + # self.assertEqual(atom_map[2], 5) + # self.assertIn(atom_map[3], [0, 1, 3, 4]) + # self.assertIn(atom_map[4], [0, 1, 3, 4]) + # self.assertIn(atom_map[5], [0, 1, 3, 4]) + # self.assertTrue(any(atom_map[r_index] in [0, 1] for r_index in [3, 4, 5])) + # self.assertIn(atom_map[6], [6, 7]) + # self.assertIn(atom_map[7], [6, 7]) + # self.assertTrue(check_atom_map(rxn)) + # + # # CH4 + OH <=> CH3 + H2O + # r_1 = ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz) + # r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) + # p_1 = ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz) + # p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[0], 0) + # self.assertIn(atom_map[1], [1, 2, 3, 5, 6]) + # self.assertIn(atom_map[2], [1, 2, 3, 5, 6]) + # self.assertIn(atom_map[3], [1, 2, 3, 5, 6]) + # self.assertIn(atom_map[4], [1, 2, 3, 5, 6]) + # self.assertEqual(atom_map[5], 4) + # self.assertIn(atom_map[6], [5, 6]) + # self.assertTrue(any(atom_map[r_index] in [5, 6] for r_index in [1, 2, 3, 4])) + # self.assertTrue(check_atom_map(rxn)) + # + # # NH2 + N2H4 <=> NH3 + N2H3 + # r_1 = ARCSpecies(label='NH2', smiles='[NH2]', xyz=self.nh2_xyz) + # r_2 = ARCSpecies(label='N2H4', smiles='NN', xyz=self.n2h4_xyz) + # p_1 = ARCSpecies(label='NH3', smiles='N', xyz=self.nh3_xyz) + # p_2 = ARCSpecies(label='N2H3', smiles='N[NH]', xyz=self.n2h3_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[0], 0) + # self.assertIn(atom_map[1], [1, 2, 3]) + # self.assertIn(atom_map[2], [1, 2, 3]) + # self.assertIn(atom_map[3], [4, 5]) + # self.assertIn(atom_map[4], [4, 5]) + # self.assertTrue(any(atom_map[r_index] in [1, 2, 3] for r_index in [5, 6, 7, 8])) + # self.assertTrue(check_atom_map(rxn)) + # + # # NH2 + N2H4 <=> N2H3 + NH3 (reversed product order compared to the above reaction) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[0], 5) + # self.assertIn(atom_map[1], [6, 7, 8]) + # self.assertIn(atom_map[2], [6, 7, 8]) + # self.assertIn(atom_map[3], [0, 1]) + # self.assertIn(atom_map[4], [0, 1]) + # self.assertTrue(any(atom_map[r_index] in [6, 7, 8] for r_index in [5, 6, 7, 8])) + # self.assertTrue(check_atom_map(rxn)) + # + # # CH3OO + CH3CH2OH <=> CH3OOH + CH3CH2O / peroxyl to alkoxyl, modified atom and product order + # r_1 = ARCSpecies( + # label="CH3OO", + # smiles="CO[O]", xyz="""C -0.41690000 0.03757000 0.00590000 + # O 0.83973000 0.69383000 -0.05239000 + # O 1.79663000 -0.33527000 -0.02406000 + # H -0.54204000 -0.62249000 -0.85805000 + # H -1.20487000 0.79501000 -0.01439000 + # H -0.50439000 -0.53527000 0.93431000""") + # r_2 = ARCSpecies(label='CH3CH2OH', smiles='CCO', xyz="""C -0.97459464 0.29181710 0.10303882 + # C 0.39565894 -0.35143697 0.10221676 + # H -1.68942501 -0.32359616 0.65926091 + # H -0.93861751 1.28685508 0.55523033 + # H -1.35943743 0.38135479 -0.91822428 + # H 0.76858330 -0.46187184 1.12485643 + # H 1.10301149 0.25256708 -0.47388355 + # O 0.30253309 -1.63748710 -0.49196889 + # H 1.19485981 -2.02360458 -0.47786539""") + # p_1 = ARCSpecies(label='CH3OOH', smiles='COO', xyz="""C -0.76039072 0.01483858 -0.00903344 + # H -1.56632337 0.61401630 -0.44251282 + # H -1.02943316 -0.30449156 1.00193709 + # O 0.16024511 1.92327904 0.86381800 + # H -0.60052507 -0.86954495 -0.63086438 + # O 0.44475333 0.76952102 0.02291303 + # H 0.30391344 2.59629139 0.17435159""") + # p_2 = ARCSpecies(label='CH3CH2O', smiles='CC[O]', xyz="""C 0.79799272 -0.01511040 0.00517437 + # H -1.13881231 -0.99286049 0.06963185 + # O 1.17260343 -0.72227959 -1.04851579 + # H -1.14162013 0.59700303 0.84092854 + # H -1.13266865 0.46233725 -0.93283228 + # C -0.74046271 0.02568566 -0.00568694 + # H 1.11374677 1.03794239 0.06905096 + # H 1.06944350 -0.38306117 1.00698657""") + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual([0,5,3],atom_map[0:3]) + # self.assertIn(tuple(atom_map[3:6]), list(permutations([1, 2, 4]))) + # self.assertEqual([12, 7], atom_map[6:8]) + # self.assertIn(tuple(atom_map[8:11]),list(permutations([8, 10, 11]))) + # self.assertIn(tuple(atom_map[11:13]),list(permutations([13, 14]))) + # self.assertEqual([9,6], atom_map[13:]) + # self.assertTrue(check_atom_map(rxn)) + # + # # C3H6O + OH <=> C3H5O + H2O + # r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) + # r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) + # p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) + # p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[:4], [0, 1, 3, 4]) + # self.assertIn(atom_map[4], [5,6, 7]) + # self.assertIn(atom_map[5], [5, 6, 7]) + # self.assertIn(atom_map[6], [5, 6, 7]) + # self.assertIn(atom_map[7], [2, 11]) + # self.assertIn(atom_map[8], [2, 11]) + # self.assertEqual(atom_map[9:], [8, 9, 10]) + # + # # C4H10O + OH <=> C4H9O + H2O + # r_1 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.c4h10o_xyz) + # r_2 = ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz) + # p_1 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) + # p_2 = ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # self.assertEqual(atom_map[:5], [0, 3, 4, 5, 6]) + # for index in [5, 6, 7]: + # self.assertIn(atom_map[index], [1, 2, 15, 16]) + # self.assertEqual(atom_map[8],7) + # for i in atom_map[9:12]: + # self.assertIn(i,[8,9,10]) + # for i in atom_map[12:14]: + # self.assertIn(i,[11,12]) + # self.assertEqual(atom_map[14],13) + # self.assertEqual(atom_map[15],14) + # self.assertIn(atom_map[16], [15, 16]) + # self.assertTrue(check_atom_map(rxn)) + # + # # C3H6O + C4H9O <=> C3H5O + C4H10O + # r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) + # r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) + # p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) + # p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.c4h10o_xyz) + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn.atom_map + # 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]) + # self.assertIn(atom_map[6], [5,6, 7]) + # self.assertIn(atom_map[7], [2, 14, 15, 16, 18, 19, 20]) + # self.assertIn(atom_map[8], [2, 14, 15, 16, 18, 19, 20]) + # self.assertIn(2, atom_map[7:9]) + # self.assertEqual(atom_map[9], 8) + # self.assertIn(atom_map[10], [9,11]) + # self.assertIn(atom_map[11], [14, 15, 16,18,19,20]) + # self.assertIn(atom_map[12], [14, 15, 16,18,19,20]) + # self.assertEqual(atom_map[13],10) + # self.assertIn(atom_map[14], [9,11]) + # self.assertEqual(atom_map[15:17], [12,13]) + # self.assertEqual(atom_map[17],17) + # self.assertIn(atom_map[18], [14, 15, 16,18,19,20]) + # self.assertIn(atom_map[19], [14, 15, 16,18,19,20]) + # self.assertIn(atom_map[20], [14, 15, 16,18,19,20]) + # self.assertIn(atom_map[21], [21,22]) + # self.assertIn(atom_map[22], [21,22]) + # self.assertEqual(atom_map[23],23) + # self.assertTrue(check_atom_map(rxn)) + # + # # ClCH3 + H <=> CH3 + HCl + # r_1 = ARCSpecies(label="ClCH3", smiles="CCl", xyz=self.ch3cl_xyz) + # r_2 = ARCSpecies(label="H", smiles="[H]", xyz=self.h_rad_xyz) + # p_1 = ARCSpecies(label="CH3", smiles="[CH3]", xyz=self.ch3_xyz_2) + # p_2 = ARCSpecies(label="HCl", smiles="[H][Cl]", xyz=self.hcl_xyz) + # rxn = ARCReaction(r_species=[r_2, r_1], p_species=[p_2, p_1]) + # atom_map = rxn.atom_map + # self.assertEqual(rxn.family.lower(),"cl_abstraction") + # self.assertEqual(atom_map[:3], [0, 1, 2]) + # for i in atom_map[3:]: + # self.assertIn(i, [3, 4, 5]) + # self.assertTrue(check_atom_map(rxn)) + # # ClCH3 + H <=> CH3 + HCl different order + # rxn_2 = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) + # atom_map = rxn_2.atom_map + # self.assertEqual(atom_map[:2], [1, 2]) + # for index in [2, 3, 4]: + # self.assertIn(atom_map[index], [3, 4, 5]) + # self.assertEqual(atom_map[-1], 0) + # self.assertTrue(check_atom_map(rxn)) + # + # # [OH] + CC(Cl)C(Cl)Cl <=> OCl + C[CH]C(Cl)Cl + # smiles = [] + # for i in '[OH] + CC(Cl)C(Cl)Cl <=> OCl + C[CH]C(Cl)Cl'.split(): + # if i != "<=>" and i != '+': + # smiles.append(i) + # + # r_1_xyz = {'symbols': ('O', 'H'), 'isotopes': (16, 1), + # 'coords': ((0.48890386738601, 0.0, 0.0), (-0.48890386738601, 0.0, 0.0))} + # + # r_2_xyz = {'symbols': ('C', 'C', 'Cl', 'C', 'Cl', 'Cl', 'H', 'H', 'H', 'H', 'H'), + # 'isotopes': (12, 12, 35, 12, 35, 35, 1, 1, 1, 1, 1), 'coords': ( + # (1.2438372893135106, 0.40661350465687324, -0.16279018264054892), + # (0.07827324125005171, -0.277154649803216, 0.5482887194488805), + # (-0.1538756923467617, 0.5009471321060629, 2.155037501334864), + # (-1.245183156820767, -0.303306879503286, -0.23533878891899096), + # (-1.1043944712471334, -1.3227416585177485, -1.7010412234762065), + # (-1.8186157680197266, 1.3177860639647956, -0.7221760707038685), + # (2.159163866798944, 0.32583527910226096, 0.4346504778666261), + # (1.056514815021544, 1.471768404816661, -0.33289291962920015), + # (1.4499964728678152, -0.05967057895051073, -1.131013164504492), + # (0.3717352549047681, -1.308596593192221, 0.7750989547682503), + # (-2.0374518517222544, -0.751480024679671, 0.37217669645466245))} + # + # p_1_xyz = {'symbols': ('O', 'Cl', 'H'), 'isotopes': (16, 35, 1), + # 'coords': ((-0.3223044372303026, 0.4343354356368888, 0.0), + # (1.2650242694442462, -0.12042710381137228, 0.0), + # (-0.9427198322139436, -0.3139083318255167, 0.0))} + # + # p_2_xyz = {'symbols': ('C', 'C', 'C', 'Cl', 'Cl', 'H', 'H', 'H', 'H', 'H'), + # 'isotopes': (12, 12, 12, 35, 35, 1, 1, 1, 1, 1), + # 'coords': ((-1.3496376883278178, -0.020445981649800302, -0.1995184115269273), + # (-0.051149096449292386, -0.3885500107837139, 0.4222976979623008), + # (1.217696701041357, 0.15947991928242372, -0.1242718714010236), + # (1.7092794464102241, 1.570982412202936, 0.8295196720275746), + # (2.474584210365428, -1.0919019396606517, -0.06869614478411318), + # (-1.6045061896547035, 1.0179450876989615, 0.03024632893682861), + # (-1.3137314500783486, -0.14754777860704252, -1.2853589013330937), + # (-2.1459595425475264, -0.6625965540242661, 0.188478021031359), + # (-0.044412318929613885, -0.9093853981117669, 1.373599947353138), + # (1.1078359281702537, 0.47202024365290884, -1.1662963382659064))} + # r_1 = ARCSpecies(label='r1', smiles=smiles[0],xyz=r_1_xyz ) + # r_2 = ARCSpecies(label='r2', smiles=smiles[1],xyz=r_2_xyz) + # p_1 = ARCSpecies(label='p1', smiles=smiles[2],xyz=p_1_xyz) + # p_2 = ARCSpecies(label='p2', smiles=smiles[3],xyz=p_2_xyz) + # + # rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn1.atom_map + # # expected: [0, 2, 3, 4, 1, 5, [6, 7], [6, 7], [8, 9, 10], [8, 9, 10], [8, 9, 10], 11, 12] + # self.assertEqual(atom_map[:6], [0,2,3,4,1,5]) + # self.assertIn(atom_map[6],[6,7]) + # self.assertIn(atom_map[7], [6, 7]) + # self.assertIn(atom_map[8], [8,9,10]) + # self.assertIn(atom_map[9], [8,9,10]) + # self.assertIn(atom_map[10], [8,9,10]) + # self.assertEqual(atom_map[11],11) + # self.assertEqual(atom_map[12], 12) + # self.assertTrue(check_atom_map(rxn)) + # + # # Br abstraction + # # OH + CH3Br <=> HOBr + CH3 + # r_1_xyz = {'symbols': ('O', 'H'), 'isotopes': (16, 1), + # 'coords': ((0.48890386738601, 0.0, 0.0), (-0.48890386738601, 0.0, 0.0))} + # + # r_2_xyz = {'symbols': ('C', 'Br', 'H', 'H', 'H'), 'isotopes': (12, 79, 1, 1, 1), 'coords': ( + # (-0.18386469024502916, -0.0018692264481234688, 0.0013619971891954718), + # (1.7508998155803106, 0.017800204658373744, -0.01296995950979447), + # (-0.5218757573028803, -0.6458197160504338, -0.8118262063895171), + # (-0.5338693855859405, 1.0212985296781085, -0.14294057406667127), + # (-0.5112899824464621, -0.3914097918379277, 0.9663747427767874))} + # + # p_1_xyz = {'symbols': ('O', 'Br', 'H'), 'isotopes': (16, 79, 1), 'coords': ( + # (-0.3691040522383542, 0.44403140947953346, 0.0), (1.3490312999095744, -0.1319682267704319, 0.0), + # (-0.9799272476712202, -0.31206318270910166, 0.0))} + # + # p_2_xyz = {'symbols': ('C', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1), 'coords': ( + # (3.3746019998564553e-09, 5.828827384106545e-09, -4.859105107686622e-09), + # (1.0669051052331406, -0.17519582095514982, 0.05416492980439295), + # (-0.6853171627400634, -0.8375353626879753, -0.028085652887100996), + # (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608))} + # + # r_1 = ARCSpecies(label='r1', smiles='[O][H]', xyz=r_1_xyz) + # r_2 = ARCSpecies(label='r2', smiles='[CH3]Br', xyz=r_2_xyz) + # p_1 = ARCSpecies(label='p1', smiles='OBr', xyz=p_1_xyz) + # p_2 = ARCSpecies(label='p2', smiles='[CH3]', xyz=p_2_xyz) + # + # rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn1.atom_map + # self.assertEqual(atom_map[:4],[0,2,3,1]) + # self.assertIn(atom_map[4], [4,5,6]) + # self.assertIn(atom_map[5], [4, 5, 6]) + # self.assertIn(atom_map[6], [4, 5, 6]) + # self.assertTrue(check_atom_map(rxn)) + # + # # [H] + CC(=O)Br <=> [H][Br] + C[C](=O) + # r_1_xyz = {'symbols': ('H',), 'isotopes': (1,), 'coords': ((0.0, 0.0, 0.0),)} + # + # r_2_xyz = {'symbols': ('C', 'C', 'O', 'Br', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 79, 1, 1, 1), + # 'coords': ((-0.7087772076387326, -0.08697184565826255, 0.08295914062572969), + # (0.7238141593293749, 0.2762480677183181, -0.14965326856248656), + # (1.1113560248255752, 1.3624373452907719, -0.554840372311578), + # (2.0636725443687616, -1.041297021241265, 0.20693447296577364), + # (-0.9844931733249197, -0.9305935329026733, -0.5546432084044857), + # (-0.8586221633621384, -0.3455305862905263, 1.134123935245044), + # (-1.3469501841979155, 0.7657075730836449, -0.16488069955797996))} + # + # p_1_xyz = {'symbols': ('C', 'C', 'O', 'H', 'H', 'H'), 'isotopes': (12, 12, 16, 1, 1, 1), + # 'coords': ((-0.4758624005470258, 0.015865899777425058, -0.11215987340300927), + # (0.9456990856850401, -0.031530842469194666, 0.2228995599390481), + # (2.0897646616994816, -0.06967555524967288, 0.492553667108967), + # (-1.08983188764878, -0.06771143046366379, 0.7892594299969324), + # (-0.7261604551815313, 0.9578749227991876, -0.6086176800339509), + # (-0.7436090040071672, -0.8048229943940851, -0.7839351036079769))} + # + # p_2_xyz = {'symbols': ('Br', 'H'), 'isotopes': (79, 1), + # 'coords': ((0.7644788559644482, 0.0, 0.0), (-0.7644788559644482, 0.0, 0.0))} + # + # r_1 = ARCSpecies(label='r1', smiles='[H]', xyz=r_1_xyz) + # r_2 = ARCSpecies(label='r2', smiles='CC(=O)Br', xyz=r_2_xyz) + # p_1 = ARCSpecies(label='p1', smiles='C[C](=O)', xyz=p_1_xyz) + # p_2 = ARCSpecies(label='p2', smiles='[Br][H]', xyz=p_2_xyz) + # + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_2, p_1]) + # atom_map=rxn.atom_map + # self.assertEqual(atom_map[:5], [1, 2, 3, 4, 0]) + # self.assertIn(tuple(atom_map[5:]), permutations([5, 6, 7])) + # self.assertTrue(check_atom_map(rxn)) + # + # # Change Order [H] + CC(=O)Br <=> C[C](=O) + [H][Br] + # r_1 = ARCSpecies(label='r1', smiles='[H]', xyz=r_1_xyz) + # r_2 = ARCSpecies(label='r2', smiles='CC(=O)Br', xyz=r_2_xyz) + # p_1 = ARCSpecies(label='p1', smiles='C[C](=O)', xyz=p_1_xyz) + # p_2 = ARCSpecies(label='p2', smiles='[H][Br]', xyz=p_2_xyz) + # + # rxn = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map=rxn.atom_map + # self.assertEqual(atom_map[:5], [7, 0, 1, 2, 6]) + # self.assertIn(tuple(atom_map[5:]), list(permutations([3, 4, 5]))) + # self.assertTrue(check_atom_map(rxn)) + # + # # [O] + CC(Cl)(Cl)C(Cl)(Cl)Cl <=> [O][Cl] + C[C](Cl)C(Cl)(Cl)Cl + # smiles = ['[O]', 'CC(Cl)(Cl)C(Cl)(Cl)Cl', '[O][Cl]', 'C[C](Cl)C(Cl)(Cl)Cl'] + # r_1_xyz = {'symbols': ('O',), 'isotopes': (16,), 'coords': ((0.0, 0.0, 0.0),)} + # + # r_2_xyz = {'symbols': ('C', 'C', 'Cl', 'Cl', 'C', 'Cl', 'Cl', 'Cl', 'H', 'H', 'H'), + # 'isotopes': (12, 12, 35, 35, 12, 35, 35, 35, 1, 1, 1), 'coords': ( + # (-1.3340513332954889, 0.2811635614535751, -0.078045907046801), + # (-0.06460593375936133, -0.5810773314093911, -0.02962891425941322), + # (-0.2609310384494481, -1.7354943987581986, 1.3623405448734305), + # (-0.06523629769352735, -1.6097818007913829, -1.5298182298699716), + # (1.2568349080206898, 0.251354210359208, 0.09596787533379413), + # (2.7373740437547514, -0.7858820942054363, 0.1510602855327231), + # (1.4729373085674606, 1.396702908938121, -1.2920641361183987), + # (1.2776463867390788, 1.2712465700052025, 1.5941477468638563), + # (-1.3327512075949484, 0.9633461541030465, -0.9346702675682734), + # (-2.235286345856216, -0.338363905821591, -0.1659562352150731), + # (-1.45193049043298, 0.886786126126846, 0.8266672374741411))} + # + # p_1_xyz = {'symbols': ('O', 'Cl'), 'isotopes': (16, 35), + # 'coords': ((0.8407400963991551, 0.0, 0.0), (-0.8407400963991551, 0.0, 0.0))} + # + # p_2_xyz = {'symbols': ('C', 'C', 'Cl', 'C', 'Cl', 'Cl', 'Cl', 'H', 'H', 'H'), + # 'isotopes': (12, 12, 35, 12, 35, 35, 35, 1, 1, 1), 'coords': ( + # (-1.3826664358998055, -0.04852445131046896, -0.016935550260331302), + # (-0.01984344739858957, 0.5351447284412386, 0.14069644461529232), + # (0.06780252918727915, 2.0178457939896477, 1.0316373428560468), + # (1.240695333262242, -0.22627953918952265, -0.15010504208991474), + # (2.5003017492701316, 0.8385176202279041, -0.8511606324628386), + # (1.8619474142609682, -0.9616513146239644, 1.3591396432655138), + # (0.9630230000989414, -1.5484613928720057, -1.3347069863893728), + # (-1.4535219021739985, -1.0095075283181074, 0.502205010423143), + # (-2.1607091682952886, 0.6031752006499635, 0.39420249485619346), + # (-1.6170290723118037, -0.20025911699469934, -1.0749727248137075))} + # + # r_1 = ARCSpecies(label='r1', smiles=smiles[0], xyz=r_1_xyz) + # r_2 = ARCSpecies(label='r2', smiles=smiles[1], xyz=r_2_xyz) + # p_1 = ARCSpecies(label='p1', smiles=smiles[2], xyz=p_1_xyz) + # p_2 = ARCSpecies(label='p2', smiles=smiles[3], xyz=p_2_xyz) + # + # rxn1 = ARCReaction(r_species=[r_1, r_2], p_species=[p_1, p_2]) + # atom_map = rxn1.atom_map + # self.assertEqual(atom_map[:3],[0,2,3]) + # self.assertIn(atom_map[3:5],[[1,4],[4,1]]) + # self.assertEqual(atom_map[5],5) + # self.assertIn(atom_map[6], [6,7,8]) + # self.assertIn(atom_map[7], [6, 7, 8]) + # self.assertIn(atom_map[8], [6, 7, 8]) + # self.assertIn(atom_map[9], [9, 10, 11]) + # self.assertIn(atom_map[10], [9, 10, 11]) + # self.assertIn(atom_map[11], [9, 10, 11]) + # self.assertTrue(check_atom_map(rxn1)) def test_map_ho2_elimination_from_peroxy_radical(self): """Test the map_ho2_elimination_from_peroxy_radical() function.""" diff --git a/arc/mapping/engine.py b/arc/mapping/engine.py index 45d97ce724..893001c6bc 100644 --- a/arc/mapping/engine.py +++ b/arc/mapping/engine.py @@ -17,9 +17,9 @@ from rmgpy.molecule import Molecule from rmgpy.species import Species -import arc.rmgdb as rmgdb 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.reaction.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 @@ -27,201 +27,268 @@ from numpy import unique if TYPE_CHECKING: - from rmgpy.data.kinetics.family import TemplateReaction - from rmgpy.data.rmg import RMGDatabase from rmgpy.molecule.molecule import Atom - from rmgpy.reaction import Reaction from arc.reaction import ARCReaction RESERVED_FINGERPRINT_KEYS = ['self', 'chirality', 'label'] -def get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction: 'ARCReaction', - rmg_reaction: 'TemplateReaction', - ) -> Tuple[Optional[Dict[str, int]], Optional[Dict[str, int]]]: +def get_atom_indices_of_labeled_atoms_in_a_reaction(arc_reaction: 'ARCReaction', + ) -> Tuple[Dict[str, int], Dict[str, int]]: """ - Get the RMG reaction atom labels and the corresponding 0-indexed atom indices - for all labeled atoms in a TemplateReaction. + Get the atom indices for all labeled atoms in an ARCReaction object instance. Args: arc_reaction (ARCReaction): An ARCReaction object instance. - rmg_reaction (TemplateReaction): A respective RMG family TemplateReaction object instance. + + Todo: + - Currently only considering the first reaction for a single atom map, + should be extended to include all reactions and make atom_map a list of ARCReaction (call it .atom_maps ?). Returns: - Tuple[Optional[Dict[str, int]], Optional[Dict[str, int]]]: - The tuple entries relate to reactants and products. - Keys are labels (e.g., '*1'), values are corresponding 0-indices atoms. - """ - if not hasattr(rmg_reaction, 'labeled_atoms') or not rmg_reaction.labeled_atoms: - return None, None - - for spc in rmg_reaction.reactants + rmg_reaction.products: - generate_resonance_structures(object_=spc, save_order=True) - - r_map, p_map = map_arc_rmg_species(arc_reaction=arc_reaction, rmg_reaction=rmg_reaction, concatenate=False) - - reactant_index_dict, product_index_dict = dict(), dict() - reactant_atoms, product_atoms = list(), list() - rmg_reactant_order = [val for _, val in sorted(r_map.items(), key=lambda item: item[0])] - rmg_product_order = [val for _, val in sorted(p_map.items(), key=lambda item: item[0])] - for i in rmg_reactant_order: - reactant_atoms.extend([atom for atom in rmg_reaction.reactants[i].atoms]) - for i in rmg_product_order: - product_atoms.extend([atom for atom in rmg_reaction.products[i].atoms]) - - for labeled_atom_dict, atom_list, index_dict in zip([rmg_reaction.labeled_atoms['reactants'], - rmg_reaction.labeled_atoms['products']], - [reactant_atoms, product_atoms], - [reactant_index_dict, product_index_dict]): - for label, atom_1 in labeled_atom_dict.items(): - for i, atom_2 in enumerate(atom_list): - if atom_1.id == atom_2.id: - index_dict[label] = i - break + Tuple[Dict[str, int], Dict[str, int]]: Keys are labels (e.g., '*1'), + values are corresponding 0-indexed atom indices + in the reactants and in the products. + """ + product_dicts = get_reaction_family_products(arc_reaction) + product_dict = product_dicts[0] + reactant_index_dict, product_index_dict = product_dict['label_map'], dict() + pairs = pair_reaction_products(arc_reaction, product_dict['products']) + maps = dict() + for arc_rxn_idx, prod_idx in pairs.items(): + maps[prod_idx] = map_two_species(arc_reaction.p_species[arc_rxn_idx], + ARCSpecies(label=f'P{prod_idx}', mol=product_dict['products'][prod_idx]), + map_type='dict', + consider_chirality=False, + ) + fam_prods_map_to_rxn_prods = dict() + for i in range(len(maps.keys())): + for key, val in maps[i].items(): + fam_prods_map_to_rxn_prods[key + sum([len(mol.atoms) for mol in product_dict['products'][:i]])] = ( + val + sum([len(spc.mol.atoms) for spc in arc_reaction.r_species[:i]])) + product_index_dict = {key: fam_prods_map_to_rxn_prods[val] for key, val in reactant_index_dict.items()} return reactant_index_dict, product_index_dict -def map_arc_rmg_species(arc_reaction: 'ARCReaction', - rmg_reaction: Union['Reaction', 'TemplateReaction'], - concatenate: bool = True, - ) -> Tuple[Dict[int, Union[List[int], int]], Dict[int, Union[List[int], int]]]: +def pair_reaction_products(reaction: 'ARCReaction', + products: List['Molecule'], + ) -> Dict[int, int]: """ - Map the species pairs in an ARC reaction to those in a respective RMG reaction - which is defined in the same direction. + Map the species pairs in an ARC reaction to those in the given product list. + This function assumes that resonance structures (mol_list) were generated for the ARCReaction object instance. Args: - arc_reaction (ARCReaction): An ARCReaction object instance. - rmg_reaction (Union[Reaction, TemplateReaction]): A respective RMG family TemplateReaction object instance. - concatenate (bool, optional): Whether to return isomorphic species as a single list (``True``, default), - or to return isomorphic species separately (``False``). + reaction (ARCReaction): An ARCReaction object instance. + products (List[ARCSpecies]): Species that correspond to the ARCReaction products that require pairing. Returns: - Tuple[Dict[int, Union[List[int], int]], Dict[int, Union[List[int], int]]]: - The first tuple entry refers to reactants, the second to products. - Keys are specie indices in the ARC reaction, - values are respective indices in the RMG reaction. - If ``concatenate`` is ``True``, values are lists of integers. Otherwise, values are integers. - """ - if rmg_reaction.is_isomerization(): - if concatenate: - return {0: [0]}, {0: [0]} - else: - return {0: 0}, {0: 0} - r_map, p_map = dict(), dict() - arc_reactants, arc_products = arc_reaction.get_reactants_and_products(arc=True) - for spc_map, rmg_species, arc_species in [(r_map, rmg_reaction.reactants, arc_reactants), - (p_map, rmg_reaction.products, arc_products)]: - for i, arc_spc in enumerate(arc_species): - for j, rmg_obj in enumerate(rmg_species): - rmg_spc = Species(molecule=[rmg_obj]) if isinstance(rmg_obj, Molecule) else rmg_obj - if not isinstance(rmg_spc, Species): - raise ValueError(f'Expected an RMG object instances of Molecule or Species, ' - f'got {rmg_obj} which is a {type(rmg_obj)}.') - generate_resonance_structures(object_=rmg_spc, save_order=True) - rmg_spc_based_on_arc_spc = Species(molecule=arc_spc.mol_list) - generate_resonance_structures(object_=rmg_spc_based_on_arc_spc, save_order=True) - if rmg_spc.is_isomorphic(rmg_spc_based_on_arc_spc, save_order=True): - if i in spc_map.keys() and concatenate: - spc_map[i].append(j) - elif concatenate: - spc_map[i] = [j] - elif i not in spc_map.keys() and j not in spc_map.values(): - spc_map[i] = j + Dict[int, int]: Keys are species indices in the ARC reaction, values are respective indices in the product list. + """ + if reaction.is_isomerization(): + return {0: 0} + product_pairs = dict() + arc_products = reaction.get_reactants_and_products(arc=True)[1] + for i, arc_product in enumerate(arc_products): + found = False + for j, prod in enumerate(products): + if j not in product_pairs.values(): + for arc_mol in arc_product.mol_list: + if arc_mol.is_isomorphic(prod): + product_pairs[i] = j + found = True break - if not r_map or not p_map: - raise ValueError(f'Could not match some of the RMG Reaction {rmg_reaction} to the ARC Reaction {arc_reaction}.') - return r_map, p_map - - -def find_equivalent_atoms_in_reactants(arc_reaction: 'ARCReaction', - backend: str = 'ARC', - ) -> List[List[int]]: - """ - Find atom indices that are equivalent in the reactants of an ARCReaction - in the sense that they represent degenerate reaction sites that are indifferentiable in 2D. - Bridges between RMG reaction templates and ARC's 3D TS structures. - Running indices in the returned structure relate to reactant_0 + reactant_1 + ... - - Args: - arc_reaction ('ARCReaction'): The ARCReaction object instance. - backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. - - Returns: - List[List[int]]: Entries are lists of 0-indices, each such list represents equivalent atoms. - """ - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction, backend=backend) - dicts = [get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(rmg_reaction=rmg_reaction, - arc_reaction=arc_reaction)[0] - for rmg_reaction in rmg_reactions] - equivalence_map = dict() - for index_dict in dicts: - for key, value in index_dict.items(): - if key in equivalence_map: - equivalence_map[key].append(value) - else: - equivalence_map[key] = [value] - equivalent_indices = list(list(set(equivalent_list)) for equivalent_list in equivalence_map.values()) - return equivalent_indices - - -def get_rmg_reactions_from_arc_reaction(arc_reaction: 'ARCReaction', - backend: str = 'ARC', - db: Optional['RMGDatabase'] = None, - ) -> Optional[List['TemplateReaction']]: - """ - A helper function for getting RMG reactions from an ARC reaction. - This function calls ``map_two_species()`` so that each species in the RMG reaction is correctly mapped - to the corresponding species in the ARC reaction. It does not attempt to map reactants to products. - - Args: - arc_reaction (ARCReaction): The 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[TemplateReaction]]: - The respective RMG TemplateReaction object instances (considering resonance structures). - """ - if arc_reaction.family is None: - rmgdb.determine_family(reaction=arc_reaction, db=db) - if arc_reaction.family is None: - return None - if not arc_reaction.family.save_order: - raise ValueError('Must have the save_order attribute of the family set to True.') - rmg_reactions = arc_reaction.family.generate_reactions(reactants=[spc.mol.copy(deep=True) for spc in arc_reaction.r_species], - products=[spc.mol.copy(deep=True) for spc in arc_reaction.p_species], - prod_resonance=True, - delete_labels=False, - relabel_atoms=False, - ) - for rmg_reaction in rmg_reactions: - r_map, p_map = map_arc_rmg_species(arc_reaction=arc_reaction, rmg_reaction=rmg_reaction, concatenate=False) - try: - ordered_rmg_reactants = [rmg_reaction.reactants[r_map[i]] for i in range(len(rmg_reaction.reactants))] - ordered_rmg_products = [rmg_reaction.products[p_map[i]] for i in range(len(rmg_reaction.products))] - except KeyError: - logger.warning(f'Got a problematic RMG rxn from ARC rxn, trying again') - continue - mapped_rmg_reactants, mapped_rmg_products = list(), list() - for ordered_rmg_mols, arc_species, mapped_mols in zip([ordered_rmg_reactants, ordered_rmg_products], - [arc_reaction.r_species, arc_reaction.p_species], - [mapped_rmg_reactants, mapped_rmg_products], - ): - for rmg_mol, arc_spc in zip(ordered_rmg_mols, arc_species): - mol = arc_spc.copy().mol - # The RMG molecule will get a random 3D conformer, don't consider chirality when mapping. - atom_map = map_two_species(mol, rmg_mol, map_type='dict', backend=backend, consider_chirality=False) - if atom_map is None: - continue - new_atoms_list = list() - for i in range(len(rmg_mol.atoms)): - rmg_mol.atoms[atom_map[i]].id = mol.atoms[i].id - new_atoms_list.append(rmg_mol.atoms[atom_map[i]]) - rmg_mol.atoms = new_atoms_list - mapped_mols.append(rmg_mol) - rmg_reaction.reactants, rmg_reaction.products = mapped_rmg_reactants, mapped_rmg_products - return rmg_reactions + if found: + break + if len(product_pairs.keys()) != len(arc_products): + raise ValueError(f'Could not match some of the products in: {reaction}\nto the given products:\n{products}') + return product_pairs + + + + + + + +# def get_rmg_reactions_from_arc_reaction(arc_reaction: 'ARCReaction', +# backend: str = 'ARC', +# ) -> Optional[List['TemplateReaction']]: +# """ +# A helper function for getting RMG reactions from an ARC reaction. +# This function calls ``map_two_species()`` so that each species in the RMG reaction is correctly mapped +# to the corresponding species in the ARC reaction. It does not attempt to map reactants to products. +# +# Args: +# arc_reaction (ARCReaction): The ARCReaction object instance. +# backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. +# +# Returns: +# Optional[List[TemplateReaction]]: +# The respective RMG TemplateReaction object instances (considering resonance structures). +# """ +# if arc_reaction.family is None: +# return None +# rmg_reactions = arc_reaction.family.generate_reactions(reactants=[spc.mol.copy(deep=True) for spc in arc_reaction.r_species], +# products=[spc.mol.copy(deep=True) for spc in arc_reaction.p_species], +# prod_resonance=True, +# delete_labels=False, +# relabel_atoms=False, +# ) +# for rmg_reaction in rmg_reactions: +# r_map, p_map = map_arc_rmg_species(arc_reaction=arc_reaction, rmg_reaction=rmg_reaction, concatenate=False) +# try: +# ordered_rmg_reactants = [rmg_reaction.reactants[r_map[i]] for i in range(len(rmg_reaction.reactants))] +# ordered_rmg_products = [rmg_reaction.products[p_map[i]] for i in range(len(rmg_reaction.products))] +# except KeyError: +# logger.warning(f'Got a problematic RMG rxn from ARC rxn, trying again') +# continue +# mapped_rmg_reactants, mapped_rmg_products = list(), list() +# for ordered_rmg_mols, arc_species, mapped_mols in zip([ordered_rmg_reactants, ordered_rmg_products], +# [arc_reaction.r_species, arc_reaction.p_species], +# [mapped_rmg_reactants, mapped_rmg_products], +# ): +# for rmg_mol, arc_spc in zip(ordered_rmg_mols, arc_species): +# mol = arc_spc.copy().mol +# # The RMG molecule will get a random 3D conformer, don't consider chirality when mapping. +# atom_map = map_two_species(mol, rmg_mol, map_type='dict', backend=backend, consider_chirality=False) +# if atom_map is None: +# continue +# new_atoms_list = list() +# for i in range(len(rmg_mol.atoms)): +# rmg_mol.atoms[atom_map[i]].id = mol.atoms[i].id +# new_atoms_list.append(rmg_mol.atoms[atom_map[i]]) +# rmg_mol.atoms = new_atoms_list +# mapped_mols.append(rmg_mol) +# rmg_reaction.reactants, rmg_reaction.products = mapped_rmg_reactants, mapped_rmg_products +# return rmg_reactions +# +# +# def get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction: 'ARCReaction', +# rmg_reaction: 'TemplateReaction', +# ) -> Tuple[Optional[Dict[str, int]], Optional[Dict[str, int]]]: +# """ +# Get the RMG reaction atom labels and the corresponding 0-indexed atom indices +# for all labeled atoms in a TemplateReaction. +# +# Args: +# arc_reaction (ARCReaction): An ARCReaction object instance. +# rmg_reaction (TemplateReaction): A respective RMG family TemplateReaction object instance. +# +# Returns: +# Tuple[Optional[Dict[str, int]], Optional[Dict[str, int]]]: +# The tuple entries relate to reactants and products. +# Keys are labels (e.g., '*1'), values are corresponding 0-indices atoms. +# """ +# if not hasattr(rmg_reaction, 'labeled_atoms') or not rmg_reaction.labeled_atoms: +# return None, None +# +# for spc in rmg_reaction.reactants + rmg_reaction.products: +# generate_resonance_structures(object_=spc, save_order=True) +# +# r_map, p_map = map_arc_rmg_species(arc_reaction=arc_reaction, rmg_reaction=rmg_reaction, concatenate=False) +# +# reactant_index_dict, product_index_dict = dict(), dict() +# reactant_atoms, product_atoms = list(), list() +# rmg_reactant_order = [val for _, val in sorted(r_map.items(), key=lambda item: item[0])] +# rmg_product_order = [val for _, val in sorted(p_map.items(), key=lambda item: item[0])] +# for i in rmg_reactant_order: +# reactant_atoms.extend([atom for atom in rmg_reaction.reactants[i].atoms]) +# for i in rmg_product_order: +# product_atoms.extend([atom for atom in rmg_reaction.products[i].atoms]) +# +# for labeled_atom_dict, atom_list, index_dict in zip([rmg_reaction.labeled_atoms['reactants'], +# rmg_reaction.labeled_atoms['products']], +# [reactant_atoms, product_atoms], +# [reactant_index_dict, product_index_dict]): +# for label, atom_1 in labeled_atom_dict.items(): +# for i, atom_2 in enumerate(atom_list): +# if atom_1.id == atom_2.id: +# index_dict[label] = i +# break +# return reactant_index_dict, product_index_dict +# +# +# def map_arc_rmg_species(arc_reaction: 'ARCReaction', +# rmg_reaction: Union['Reaction', 'TemplateReaction'], +# concatenate: bool = True, +# ) -> Tuple[Dict[int, Union[List[int], int]], Dict[int, Union[List[int], int]]]: +# """ +# Map the species pairs in an ARC reaction to those in a respective RMG reaction +# which is defined in the same direction. +# +# Args: +# arc_reaction (ARCReaction): An ARCReaction object instance. +# rmg_reaction (Union[Reaction, TemplateReaction]): A respective RMG family TemplateReaction object instance. +# concatenate (bool, optional): Whether to return isomorphic species as a single list (``True``, default), +# or to return isomorphic species separately (``False``). +# +# Returns: +# Tuple[Dict[int, Union[List[int], int]], Dict[int, Union[List[int], int]]]: +# The first tuple entry refers to reactants, the second to products. +# Keys are specied indices in the ARC reaction, +# values are respective indices in the RMG reaction. +# If ``concatenate`` is ``True``, values are lists of integers. Otherwise, values are integers. +# """ +# if rmg_reaction.is_isomerization(): +# if concatenate: +# return {0: [0]}, {0: [0]} +# else: +# return {0: 0}, {0: 0} +# r_map, p_map = dict(), dict() +# arc_reactants, arc_products = arc_reaction.get_reactants_and_products(arc=True) +# for spc_map, rmg_species, arc_species in [(r_map, rmg_reaction.reactants, arc_reactants), +# (p_map, rmg_reaction.products, arc_products)]: +# for i, arc_spc in enumerate(arc_species): +# for j, rmg_obj in enumerate(rmg_species): +# rmg_spc = Species(molecule=[rmg_obj]) if isinstance(rmg_obj, Molecule) else rmg_obj +# if not isinstance(rmg_spc, Species): +# raise ValueError(f'Expected an RMG object instances of Molecule or Species, ' +# f'got {rmg_obj} which is a {type(rmg_obj)}.') +# generate_resonance_structures(object_=rmg_spc, save_order=True) +# rmg_spc_based_on_arc_spc = Species(molecule=arc_spc.mol_list) +# generate_resonance_structures(object_=rmg_spc_based_on_arc_spc, save_order=True) +# if rmg_spc.is_isomorphic(rmg_spc_based_on_arc_spc, save_order=True): +# if i in spc_map.keys() and concatenate: +# spc_map[i].append(j) +# elif concatenate: +# spc_map[i] = [j] +# elif i not in spc_map.keys() and j not in spc_map.values(): +# spc_map[i] = j +# break +# if not r_map or not p_map: +# raise ValueError(f'Could not match some of the RMG Reaction {rmg_reaction} to the ARC Reaction {arc_reaction}.') +# return r_map, p_map +# +# +# def find_equivalent_atoms_in_reactants(arc_reaction: 'ARCReaction', # ? +# backend: str = 'ARC', +# ) -> List[List[int]]: +# """ +# Find atom indices that are equivalent in the reactants of an ARCReaction +# in the sense that they represent degenerate reaction sites that are indifferentiable in 2D. +# Bridges between RMG reaction templates and ARC's 3D TS structures. +# Running indices in the returned structure relate to reactant_0 + reactant_1 + ... +# +# Args: +# arc_reaction ('ARCReaction'): The ARCReaction object instance. +# backend (str, optional): Whether to use ``'QCElemental'`` or ``ARC``'s method as the backend. +# +# Returns: +# List[List[int]]: Entries are lists of 0-indices, each such list represents equivalent atoms. +# """ +# rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction, backend=backend) +# dicts = [get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(rmg_reaction=rmg_reaction, +# arc_reaction=arc_reaction)[0] +# for rmg_reaction in rmg_reactions] +# equivalence_map = dict() +# for index_dict in dicts: +# for key, value in index_dict.items(): +# if key in equivalence_map: +# equivalence_map[key].append(value) +# else: +# equivalence_map[key] = [value] +# equivalent_indices = list(list(set(equivalent_list)) for equivalent_list in equivalence_map.values()) +# return equivalent_indices def map_two_species(spc_1: Union[ARCSpecies, Species, Molecule], @@ -1031,8 +1098,7 @@ 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: + for action in ReactionFamily(label=rxn.family).actions: if action[0].lower() == "CHANGE_BOND".lower(): indicies = r_label_dict[action[1]],r_label_dict[action[3]] for r_cut in r_cuts: @@ -1064,26 +1130,51 @@ def make_bond_changes(rxn: 'ARCReaction', r_cut.mol.update() +def get_label_dict(rxn: 'ARCReaction') -> Optional[Dict[str, int]]: + """ + A function that returns the labels of the products. + + Args: + rxn: ARCReaction object + + Returns: + Optional[Dict[str, int]]: The labels of the products. + """ + products = get_reaction_family_products(rxn=rxn) + if len(products): + return products[0]['label_map'] + return None + + def assign_labels_to_products(rxn: 'ARCReaction', - p_label_dict: dict): + products: List[Molecule], + ): """ Add the indices to the reactants and products. Args: - rxn: ARCReaction object to be mapped - p_label_dict: the labels of the products - Consider changing in rmgpy. + rxn ('ARCReaction'): The reaction to be mapped. + products (List[Molecule]): The products generated from the RMG family with the same atom order as the reactants. Returns: Adding labels to the atoms of the reactants and products, to be identified later. """ - + label_dict = get_label_dict(rxn) + atom_index = 0 + for r in rxn.r_species: + for atom in r.mol.atoms: + if atom_index in label_dict.values(): + atom.label = key_by_val(label_dict, atom_index) + atom_index += 1 + product_pairs = pair_reaction_products(reaction=rxn, products=products) atom_index = 0 - for product in rxn.p_species: - for atom in product.mol.atoms: - if atom_index in p_label_dict.values() and (atom.label is str or atom.label is None): - atom.label = key_by_val(p_label_dict,atom_index) - atom_index+=1 + for product_index in range(len(products)): + rxn_product, fam_product = rxn.p_species[product_index], products[product_pairs[product_index]] + atom_map = map_two_species(spc_1=rxn_product, spc_2=fam_product, map_type='list') + for i, atom in enumerate(fam_product.atoms): + if atom_index in label_dict.values(): + rxn_product.mol.atoms[atom_map[i]].label = key_by_val(label_dict, atom_index) + atom_index += 1 def update_xyz(spcs: List[ARCSpecies]) -> List[ARCSpecies]: @@ -1095,7 +1186,7 @@ def update_xyz(spcs: List[ARCSpecies]) -> List[ARCSpecies]: spcs: the scission products that needs to be updated Returns: - new: A newely generated copies of the ARCSpecies, with updated xyz + list: A newly generated copies of the ARCSpecies, with updated xyz. """ new = list() for spc in spcs: @@ -1125,7 +1216,7 @@ def r_cut_p_cut_isomorphic(reactant, product): def pairing_reactants_and_products_for_mapping(r_cuts: List[ARCSpecies], p_cuts: List[ARCSpecies] - )-> List[Tuple[ARCSpecies,ARCSpecies]]: + ) -> List[Tuple[ARCSpecies,ARCSpecies]]: """ A function for matching reactants and products in scissored products. @@ -1134,21 +1225,21 @@ def pairing_reactants_and_products_for_mapping(r_cuts: List[ARCSpecies], p_cuts: A list of the scissored species in the reactants Returns: - a list of paired reactant and products, to be sent to map_two_species. + list: Paired reactant and products, to be sent to map_two_species. """ pairs = [] for reactant_cut in r_cuts: for product_cut in p_cuts: if r_cut_p_cut_isomorphic(reactant_cut, product_cut): pairs.append((reactant_cut, product_cut)) - p_cuts.remove(product_cut) # Just in case there are two of the same species in the list, matching them by order. + p_cuts.remove(product_cut) # Just in case there are two of the same species in the list, matching them by order. break return pairs def map_pairs(pairs): """ - A function that maps the mached species together + A function that maps the matched species together Args: pairs: A list of the pairs of reactants and species @@ -1156,11 +1247,9 @@ def map_pairs(pairs): Returns: A list of the mapped species """ - maps = list() for pair in pairs: maps.append(map_two_species(pair[0], pair[1])) - return maps @@ -1171,11 +1260,11 @@ def label_species_atoms(spcs): Args: spcs: ARCSpecies object to be labeled. """ - index=0 + index = 0 for spc in spcs: for atom in spc.mol.atoms: atom.label = str(index) - index+=1 + index += 1 def glue_maps(maps, pairs_of_reactant_and_products): @@ -1183,11 +1272,11 @@ def glue_maps(maps, pairs_of_reactant_and_products): a function that joins together the maps from the parts of the reaction. Args: - rxn: ARCReaction that requires atom mapping maps: The list of all maps of the isomorphic cuts. + pairs_of_reactant_and_products: The pairs of the reactants and products. Returns: - an Atom Map of the compleate reaction. + list: An Atom Map of the complete reaction. """ am_dict = dict() for _map, pair in zip(maps, pairs_of_reactant_and_products): @@ -1229,18 +1318,21 @@ def determine_bdes_on_spc_based_on_atom_labels(spc: "ARCSpecies", bde: Tuple[int return False -def cut_species_based_on_atom_indices(species: List["ARCSpecies"], bdes: List[Tuple[int, int]]) -> Optional[List["ARCSpecies"]]: +def cut_species_based_on_atom_indices(species: List["ARCSpecies"], + bdes: List[Tuple[int, int]], + ) -> Optional[List["ARCSpecies"]]: """ A function for scissoring species based on their atom indices. + Args: species (List[ARCSpecies]): The species list that requires scission. bdes (List[Tuple[int, int]]): A list of the atoms between which the bond should be scissored. The atoms are described using the atom labels, and not the actuall atom positions. + Returns: Optional[List["ARCSpecies"]]: The species list input after the scission. """ if not bdes: return species - for bde in bdes: for index, spc in enumerate(species): if determine_bdes_on_spc_based_on_atom_labels(spc, bde): @@ -1261,7 +1353,6 @@ def cut_species_based_on_atom_indices(species: List["ARCSpecies"], bdes: List[Tu except SpeciesError: return None break - return species @@ -1280,18 +1371,34 @@ def copy_species_list_for_mapping(species: List["ARCSpecies"]) -> List["ARCSpeci return copies -def find_all_bdes(rxn: "ARCReaction", label_dict: dict, is_reactants: bool) -> List[Tuple[int, int]]: +def find_all_bdes(rxn: "ARCReaction", + is_reactants: bool, + products: Optional[List["Molecule"]] = None, + ) -> List[Tuple[int, int]]: """ A function for finding all the broken(/formed) bonds during a chemical reaction, based on the atom indices. + Args: rxn (ARCReaction): The reaction in question. - label_dict (dict): A dictionary of the atom indices to the atom labels. - is_reactants (bool): Whether or not the species list represents reactants or products. + is_reactants (bool): Whether the species list represents reactants or products. + products (List[Molecule], optional): The products generated from the RMG family with the same atom order + as the reactants. If given, the BDE values will be mapped from them + to the reaction products. + Returns: - List[Tuple[int, int]]: A list of tuples of the form (atom_index1, atom_index2) for each broken bond. Note that these represent the atom indicies to be cut, and not final BDEs. + List[Tuple[int, int]]: A list of tuples of the form (atom_index1, atom_index2) for each broken bond. + Note that these represent the atom indices to be cut, and not final BDEs. """ + label_dict = get_label_dict(rxn) + if not is_reactants: + product_pairs = pair_reaction_products(reaction=rxn, products=products) if products is not None else None + print(label_dict) bdes = list() - for action in rxn.family.forward_recipe.actions: - if action[0].lower() == ("break_bond" if is_reactants else "form_bond"): - bdes.append((label_dict[action[1]] + 1, label_dict[action[3]] + 1)) + if rxn.family is not None: + for action in ReactionFamily(rxn.family).actions: + print(action) + if (action[0].lower() == "break_bond" and is_reactants + or action[0].lower() == "form_bond" and not is_reactants): + print(f'appending {action[1]} and {action[3]}: {(label_dict[action[1]] + 1, label_dict[action[3]] + 1)}') + bdes.append((label_dict[action[1]] + 1, label_dict[action[3]] + 1)) return bdes diff --git a/arc/mapping/engine_test.py b/arc/mapping/engine_test.py index eaf6cd2368..f56eb4fd3c 100644 --- a/arc/mapping/engine_test.py +++ b/arc/mapping/engine_test.py @@ -7,15 +7,11 @@ import unittest import numpy as np -from rmgpy.data.rmg import RMGDatabase from random import shuffle import itertools -from arc.common import _check_r_n_p_symbols_between_rmg_and_arc_rxns -from arc.species import ARCSpecies from arc.mapping.engine import * from arc.reaction import ARCReaction -from arc.rmgdb import load_families_only, determine_family from arc.species.vectors import calculate_dihedral_angle @@ -35,68 +31,28 @@ def setUpClass(cls): """ A method that is run before all unit tests in this class. """ - cls.db = RMGDatabase() - load_families_only(cls.db, "all") cls.maxDiff = None - smiles = ['CC(C)F', '[CH3]', 'C[CH](C)', 'CF'] - - r_1_1_xyz = {'symbols': ('C', 'C', 'C', 'F', 'H', 'H', 'H', 'H', 'H', 'H', 'H'), - 'isotopes': (12, 12, 12, 19, 1, 1, 1, 1, 1, 1, 1), - 'coords': ((1.2509680857915237, 0.00832885083067477, -0.28594855682006387), - (-0.08450322338173592, -0.5786110309038947, 0.12835305965368538), - (-1.196883483105121, 0.4516770584363101, 0.10106807955582568), - (0.03212452836861426, -1.0465351442062332, 1.402047416169314), - (1.2170230403876368, 0.39373449465586885, -1.309310880313081), - (1.5446944155971303, 0.8206316657310906, 0.38700047363833845), - (2.0327466889922805, -0.7555292157466509, -0.22527487012253536), - (-0.3397419937928473, -1.4280299782557704, -0.5129583662636836), - (-0.9791793765226446, 1.2777482351478369, 0.786037216866474), - (-1.340583396165929, 0.8569620299504027, -0.9049411765144166), - (-2.1366652861689137, -0.00037696563964776297, 0.43392760415012316))} - - r_2_1_xyz = {'symbols': ('C', 'H', 'H', 'H'), - 'isotopes': (12, 1, 1, 1), - 'coords' : ((3.3746019998564553e-09, 5.828827384106545e-09, -4.859105107686622e-09), - (1.0669051052331406, -0.17519582095514982, 0.05416492980439295), - (-0.6853171627400634, -0.8375353626879753, -0.028085652887100996), - (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608))} - - p_1_1_xyz = {'symbols': ('C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H'), - 'isotopes': (12, 12, 12, 1, 1, 1, 1, 1, 1, 1), - 'coords': ((-1.288730238258946, 0.06292843803165035, 0.10889818910854648), - (0.01096160773224897, -0.45756396262445836, -0.3934214957819532), - (1.2841030977199492, 0.11324607936811129, 0.12206176848573647), - (-1.4984446521053447, 1.0458196461796345, -0.3223873567509909), - (-1.2824724918369017, 0.14649429503996203, 1.1995362776757934), - (-2.098384694966955, -0.616646552269074, -0.17318515188247927), - (0.027360233461550892, -1.0601383387124987, -1.2952225290380646), - (2.122551165381095, -0.534098313164123, -0.15158596254231563), - (1.2634262459696732, 0.19628891975881263, 1.2125616721427255), - (1.4596297269035956, 1.1036697883919826, -0.307255411416999))} - - p_2_1_xyz = {'symbols': ('C', 'F', 'H', 'H', 'H'), - 'isotopes': (12, 19, 1, 1, 1), - 'coords': ((-0.060384822736851786, 0.004838867136375763, -0.004814368798794687), - (1.2877092002693546, -0.10318918150563985, 0.10266661058725791), - (-0.2965861926821434, 0.9189121874074381, -0.5532990701789506), - (-0.44047773762823295, -0.8660709320146035, -0.5425894744224189), - (-0.49026044722212864, 0.04550905897643097, 0.9980363028129072))} - cls.r_1 = ARCSpecies(label='r1', smiles=smiles[0],xyz=r_1_1_xyz ) - cls.r_2 = ARCSpecies(label='r2', smiles=smiles[1],xyz=r_2_1_xyz) - cls.p_1 = ARCSpecies(label='p1', smiles=smiles[2],xyz=p_1_1_xyz) - cls.p_2 = ARCSpecies(label='p2', smiles=smiles[3],xyz=p_2_1_xyz) - - cls.r_1_1 = ARCSpecies(label='r1', smiles=smiles[0],xyz=r_1_1_xyz ) - cls.r_2_2 = ARCSpecies(label='r2', smiles=smiles[1],xyz=r_2_1_xyz) - cls.p_1_1 = ARCSpecies(label='p1', smiles=smiles[2],xyz=p_1_1_xyz) - cls.p_2_2 = ARCSpecies(label='p2', smiles=smiles[3],xyz=p_2_1_xyz) - - cls.rxn_1 = ARCReaction(r_species=[cls.r_1, cls.r_2], p_species=[cls.p_1, cls.p_2]) - cls.rxn_1.determine_family(cls.db) - cls.rmg_reactions_rxn_1 = get_rmg_reactions_from_arc_reaction(arc_reaction=cls.rxn_1, backend="ARC") - cls.r_label_dict_rxn_1, cls.p_label_dict_rxn_1 = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=cls.rxn_1, - rmg_reaction=cls.rmg_reactions_rxn_1[0]) + smiles = ['[C]#CF', 'C#N', 'C#CF', '[C]#N'] + r_1_xyz = {'symbols': ('C', 'C', 'F'), 'isotopes': (12, 12, 19), + 'coords': ((-1.4743294980458657, -0.0015578181856664775, 0.0), + (0.06566963887880228, 6.93883492559459e-05, 0.0), + (1.4086598591670647, 0.0014884298364102333, 0.0))} + r_2_xyz = {'symbols': ('C', 'N', 'H'), 'isotopes': (12, 14, 1), + 'coords': ((-0.03166549265366707, 0.00027268185540903927, 0.0), + (1.1282915045547195, -0.009715961185409256, 0.0), + (-1.096626011901052, 0.009443279330001001, 0.0))} + p_1_xyz = {'symbols': ('C', 'C', 'F', 'H'), 'isotopes': (12, 12, 19, 1), + 'coords': ((-0.6556450958206512, 0.008288859693035082, 0.041287880045658594), + (0.5194805538847015, -0.006567413538153382, 0.28295662316778103), + (1.8344056155700479, -0.023191047733068296, 0.5533756165580456), + (-1.6982410736345843, 0.0214696015778594, -0.17312568563836497))} + p_2_xyz = {'symbols': ('C', 'N'), 'isotopes': (12, 14), 'coords': ((0.0, 0.0, 0.735), (0.0, 0.0, -0.735))} + + cls.r_1 = ARCSpecies(label='r1', smiles=smiles[0], xyz=r_1_xyz) + cls.r_2 = ARCSpecies(label='r2', smiles=smiles[1], xyz=r_2_xyz) + cls.p_1 = ARCSpecies(label='p1', smiles=smiles[2], xyz=p_1_xyz) + cls.p_2 = ARCSpecies(label='p2', smiles=smiles[3], xyz=p_2_xyz) cls.spc1 = ARCSpecies(label="Test_is_isomorphic_1",smiles="C(CO)CC") cls.spc2 = ARCSpecies(label="Test_is_isomorphic_2",smiles="OCCCC") @@ -134,15 +90,12 @@ def setUpClass(cls): (0.7709290243761263, -0.38422053705817527, 1.4054470816682596), (1.337910696892115, -1.3171321272490044, 0.001256204546378134), (2.595292874972368, 0.5254618254772234, 0.4066018700054956))}) - - cls.ch4_xyz = {'symbols': ('C', 'H', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1, 1), 'coords': ((-5.45906343962835e-10, 4.233517924761169e-10, 2.9505240956083194e-10), (-0.6505520089868748, -0.7742801979689132, -0.4125187934483119), (-0.34927557824779626, 0.9815958255612931, -0.3276823191685369), (-0.022337921721882443, -0.04887374527620588, 1.0908766524267022), (1.0221655095024578, -0.15844188273952128, -0.350675540104908))} - cls.ch3cl_xyz = { "symbols": ("Cl", "C", "H", "H", "H"), "isotopes": (35, 12, 1, 1, 1), @@ -151,10 +104,7 @@ def setUpClass(cls): (-0.13902319255720844, -0.005942522549225256, 0.0025040105508790877), (-0.48751519000626353, -0.5782116958350602, -0.8600292104425608), (-0.45894137315516464, -0.4942789247294056, 0.9255869621295756), - (-0.540565317573775, 1.0089280316586664, -0.038774361328407066), - ), - } - + (-0.540565317573775, 1.0089280316586664, -0.038774361328407066))} cls.ch3_xyz_2 = { "symbols": ("C", "H", "H", "H"), "isotopes": (12, 1, 1, 1), @@ -162,18 +112,12 @@ def setUpClass(cls): (3.3746019998564553e-09, 5.828827384106545e-09, -4.859105107686622e-09), (1.0669051052331406, -0.17519582095514982, 0.05416492980439295), (-0.6853171627400634, -0.8375353626879753, -0.028085652887100996), - (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608), - ), - } - + (-0.3815879458676787, 1.0127311778142964, -0.026079272058187608))} cls.h_rad_xyz = {"symbols": ("H",), "isotopes": (1,), "coords": ((0, 0, 0),)} - cls.hcl_xyz = { "symbols": ("H", "Cl"), "isotopes": (1, 35), - "coords": ((0.6878248644303301, 0.0, 0.0), (-0.6878248644303301, 0.0, 0.0)), - } - + "coords": ((0.6878248644303301, 0.0, 0.0), (-0.6878248644303301, 0.0, 0.0))} cls.ch4_xyz = { "symbols": ("C", "H", "H", "H", "H"), "isotopes": (12, 1, 1, 1, 1), @@ -182,10 +126,7 @@ def setUpClass(cls): (-0.6505520089868748, -0.7742801979689132, -0.4125187934483119), (-0.34927557824779626, 0.9815958255612931, -0.3276823191685369), (-0.022337921721882443, -0.04887374527620588, 1.0908766524267022), - (1.0221655095024578, -0.15844188273952128, -0.350675540104908), - ), - } - + (1.0221655095024578, -0.15844188273952128, -0.350675540104908))} cls.ch4_xyz_diff_order = """H -0.65055201 -0.77428020 -0.41251879 H -0.34927558 0.98159583 -0.32768232 C -0.00000000 0.00000000 0.00000000 @@ -541,16 +482,40 @@ def setUpClass(cls): H -1.04996634 -0.37234114 0.91874740 H 1.36260637 0.37153887 -0.86221771""" + def test_get_atom_indices_of_labeled_atoms_in_a_reaction(self): + """Test getting atom indices of labeled atoms in a reaction""" + atom_indices = get_atom_indices_of_labeled_atoms_in_a_reaction(self.arc_reaction_1) + self.assertEqual(atom_indices, ({'*1': 0, '*2': 1, '*3': 5}, {'*1': 0, '*2': 1, '*3': 7})) + + atom_indices = get_atom_indices_of_labeled_atoms_in_a_reaction(self.arc_reaction_2) + print(atom_indices) + self.assertEqual(atom_indices, ({'*1': 0, '*2': 3, '*3': 11}, {'*1': 0, '*2': 3, '*3': 16})) + + def test_pair_reaction_products(self): + """Test pairing reaction and products""" + products = [Molecule(smiles='[CH3]'), Molecule(smiles='O')] + product_pairs = pair_reaction_products(self.arc_reaction_1, products) + self.assertEqual(product_pairs, {0: 0, 1: 1}) + + product_pairs = pair_reaction_products(self.arc_reaction_1, [products[1], products[0]]) + self.assertEqual(product_pairs, {0: 1, 1: 0}) + + rxn_1 = ARCReaction(label='H2NN(T) + N2H4 <=> N2H3 + N2H3', + r_species=[ARCSpecies(label='H2NN(T)', smiles='N[N]'), + ARCSpecies(label='N2H4', smiles='NN')], + p_species=[ARCSpecies(label='N2H3', smiles='N[NH]')]) + prod_mol = Molecule(smiles='N[NH]') + product_pairs = pair_reaction_products(rxn_1, [prod_mol, prod_mol]) + self.assertEqual(product_pairs, {0: 0, 1: 1}) def test_assign_labels_to_products(self): + """Test assigning labels to products based on the atom map of the reaction""" rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - assign_labels_to_products(rxn_1_test, self.p_label_dict_rxn_1) - index = 0 - for product in rxn_1_test.p_species: - for atom in product.mol.atoms: - if not isinstance(atom.label, str) or atom.label != "": - self.assertEqual(self.p_label_dict_rxn_1[atom.label], index) - index+=1 + assign_labels_to_products(rxn_1_test, rxn_1_test.get_family_products()) + self.assertEqual([atom.label for atom in rxn_1_test.r_species[0].mol.atoms], ['*3', '', '']) + self.assertEqual([atom.label for atom in rxn_1_test.r_species[1].mol.atoms], ['*1', '', '*2']) + self.assertEqual([atom.label for atom in rxn_1_test.p_species[0].mol.atoms], ['*3', '', '*1', '']) + self.assertEqual([atom.label for atom in rxn_1_test.p_species[1].mol.atoms], ['', '*2']) def test_inc_vals(self): """Test creating an atom map via map_two_species() and incrementing all values""" @@ -561,49 +526,44 @@ def test_inc_vals(self): def test_label_species_atoms(self): rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) - assign_labels_to_products(rxn_1_test, self.p_label_dict_rxn_1) - + assign_labels_to_products(rxn_1_test) reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) - label_species_atoms(reactants) label_species_atoms(products) - index = 0 for reactant in reactants: for atom in reactant.mol.atoms: self.assertEqual(atom.label,str(index)) - index +=1 - + index += 1 index = 0 for product in products: for atom in product.mol.atoms: - self.assertEqual(atom.label,str(index)) - index +=1 + self.assertEqual(atom.label, str(index)) + index += 1 def test_cut_species_based_on_atom_indices(self): """test the cut_species_for_mapping function""" - rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) - reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) + rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2], + rmg_family_set=['H_Abstraction']) + 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_bdes(rxn_1_test, self.r_label_dict_rxn_1, True), find_all_bdes(rxn_1_test, self.p_label_dict_rxn_1, False) - + r_bdes, p_bdes = find_all_bdes(rxn_1_test, True), find_all_bdes(rxn_1_test, False) + print(r_bdes, p_bdes) r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) p_cuts = cut_species_based_on_atom_indices(products, p_bdes) + print([a.mol for a in r_cuts], [a.mol for a in p_cuts]) + self.assertIn('[C]#CF', [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) + self.assertIn('[C]#N', [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('[C]#CF', [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) + self.assertIn('[C]#N', [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) + self.assertIn('[H]', [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) - self.assertIn("C[CH]C", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) - self.assertIn("[F]", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) - self.assertIn("[CH3]", [r_cut.mol.copy(deep=True).smiles for r_cut in r_cuts]) - self.assertIn("C[CH]C", [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) - self.assertIn("[F]", [p_cut.mol.copy(deep=True).smiles for p_cut in p_cuts]) - self.assertIn("[CH3]", [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", bdes=[(1, 2), (2, 3)]) for i, a in enumerate(spc.mol.atoms): - a.label=str(i) + a.label = str(i) cuts = cut_species_based_on_atom_indices([spc], [(1, 2), (2, 3)]) self.assertEqual(len(cuts), 3) for cut in cuts: @@ -613,24 +573,21 @@ 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], [(1, 2)]) 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')] label_species_atoms(spcs) cuts = cut_species_based_on_atom_indices(spcs, [(6, 5), (4, 2), (3, 7)]) self.assertEqual(len(cuts), 4) def test_r_cut_p_cut_isomorphic(self): rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_bdes(rxn_1_test, self.r_label_dict_rxn_1, True), find_all_bdes(rxn_1_test, self.p_label_dict_rxn_1, False) - + r_bdes, p_bdes = find_all_bdes(rxn_1_test, True), find_all_bdes(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) @@ -647,12 +604,9 @@ def test_r_cut_p_cut_isomorphic(self): def test_pairing_reactants_and_products_for_mapping(self): smiles = ["[F]", "C[CH]C", "[CH3]"] rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_bdes(rxn_1_test, self.r_label_dict_rxn_1, True), find_all_bdes(rxn_1_test, self.p_label_dict_rxn_1, False) - + r_bdes, p_bdes = find_all_bdes(rxn_1_test, True), find_all_bdes(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) @@ -663,17 +617,11 @@ def test_pairing_reactants_and_products_for_mapping(self): self.assertIn(str(pair[0].mol.copy(deep=True).to_smiles()), smiles) smiles.remove(pair[0].mol.copy(deep=True).to_smiles()) - def test_map_pairs(self): + def test_map_pairs(self): # H abs reverse 1st rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) - 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) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_bdes(rxn_1_test, r_label_dict, True), find_all_bdes(rxn_1_test, p_label_dict, False) - + r_bdes, p_bdes = find_all_bdes(rxn_1_test, True), find_all_bdes(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) @@ -695,16 +643,10 @@ def test_map_pairs(self): def test_glue_maps(self): rxn_1_test = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn_1_test.determine_family(self.db) - 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) + assign_labels_to_products(rxn_1_test) reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species) label_species_atoms(reactants), label_species_atoms(products) - - r_bdes, p_bdes = find_all_bdes(rxn_1_test, r_label_dict, True), find_all_bdes(rxn_1_test, p_label_dict, False) - + r_bdes, p_bdes = find_all_bdes(rxn_1_test, True), find_all_bdes(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) @@ -716,213 +658,205 @@ def test_glue_maps(self): for index, value in enumerate(atom_map): self.assertEqual(atoms_r[index].symbol, atoms_p[value].symbol) - def test_get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(self): - """Test the get_atom_indices_of_labeled_atoms_in_an_rmg_reaction() function.""" - determine_family(self.arc_reaction_1, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_1, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_1, - rmg_reaction=rmg_reactions[0]) - self.assertEqual(r_dict['*1'], 0) - self.assertIn(r_dict['*2'], [1, 2, 3, 4]) - self.assertEqual(r_dict['*3'], 5) - self.assertEqual(p_dict['*1'], 0) - self.assertIn(p_dict['*2'], [5, 6]) - self.assertEqual(p_dict['*3'], 4) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_1, rmg_reactions)) - - determine_family(self.arc_reaction_2, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_2, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_2, - rmg_reaction=rmg_reactions[0]) - self.assertIn(r_dict['*1'], [0, 2]) - self.assertIn(r_dict['*2'], [3, 4, 5, 8, 9, 10]) - self.assertEqual(r_dict['*3'], 11) - self.assertEqual(p_dict['*1'], 0) - self.assertIn(p_dict['*2'], [11, 12, 13]) - self.assertEqual(p_dict['*3'], 10) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_2, rmg_reactions)) - - determine_family(self.arc_reaction_4, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_4, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_4, - rmg_reaction=rmg_reactions[0]) - self.assertEqual(r_dict['*1'], 0) - self.assertEqual(r_dict['*2'], 2) - self.assertIn(r_dict['*3'], [7, 8]) - self.assertEqual(p_dict['*1'], 0) - self.assertEqual(p_dict['*2'], 2) - self.assertIn(p_dict['*3'], [3, 4, 5]) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_4, rmg_reactions)) - - determine_family(self.rxn_2a, db=self.db) - for atom, symbol in zip(self.rxn_2a.r_species[0].mol.atoms, ['C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H']): - self.assertEqual(atom.symbol, symbol) - self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[0].radical_electrons, 0) - self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[1].radical_electrons, 1) - self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[2].radical_electrons, 0) - self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[0].radical_electrons, 0) - self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[1].radical_electrons, 0) - self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[2].radical_electrons, 1) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.rxn_2a, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.rxn_2a, - rmg_reaction=rmg_reactions[0]) - self.assertEqual(r_dict['*1'], 1) - self.assertIn(r_dict['*2'], [0, 2]) - self.assertIn(r_dict['*3'], [4, 5, 6, 7, 8, 9]) - self.assertEqual(p_dict['*1'], 1) - self.assertEqual(p_dict['*2'], 2) - self.assertIn(p_dict['*3'], [3, 6]) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.rxn_2a, rmg_reactions)) - - determine_family(self.rxn_2b, db=self.db) - for atom, symbol in zip(self.rxn_2b.r_species[0].mol.atoms, ['C', 'C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H']): - self.assertEqual(atom.symbol, symbol) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.rxn_2b, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.rxn_2b, - rmg_reaction=rmg_reactions[0]) - self.assertEqual(r_dict['*1'], 1) - self.assertIn(r_dict['*2'], [0, 6]) - self.assertIn(r_dict['*3'], [3, 4, 5, 7, 8, 9]) - self.assertEqual(p_dict['*1'], 1) - self.assertEqual(p_dict['*2'], 2) - self.assertIn(p_dict['*3'], [3, 6]) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.rxn_2b, rmg_reactions)) - - # C3H6O + C4H9O <=> C3H5O + C4H10O - r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O') - r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO') - p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O') - p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO') - rxn_1 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], - r_species=[r_1, r_2], p_species=[p_1, p_2]) - determine_family(rxn_1, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_1, db=self.db) - for rmg_reaction in rmg_reactions: - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_1, - rmg_reaction=rmg_reaction) - for d in [r_dict, p_dict]: - self.assertEqual(len(list(d.keys())), 3) - keys = list(d.keys()) - for label in ['*1', '*2', '*3']: - self.assertIn(label, keys) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_1, rmg_reactions)) - - p_1 = ARCSpecies(label='C3H5O', smiles='CC=C[O]') # Use a wrong resonance structure and repeat the above. - rxn_2 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], - r_species=[r_1, r_2], p_species=[p_1, p_2]) - determine_family(rxn_2, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_2, db=self.db) - for rmg_reaction in rmg_reactions: - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_2, - rmg_reaction=rmg_reaction) - for d in [r_dict, p_dict]: - self.assertEqual(len(list(d.keys())), 3) - keys = list(d.keys()) - for label in ['*1', '*2', '*3']: - self.assertIn(label, keys) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_2, rmg_reactions)) - - # C3H6O + C4H9O <=> C3H5O + C4H10O - r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) - r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) - p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) - p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.c4h10o_xyz) - rxn_3 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], - r_species=[r_1, r_2], p_species=[p_1, p_2]) - determine_family(rxn_3, db=self.db) - rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_3, db=self.db) - r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_3, - rmg_reaction=rmg_reactions[0]) - self.assertEqual(r_dict, {'*3': 10, '*1': 1, '*2': 7}) - self.assertEqual(p_dict, {'*1': 1, '*3': 9, '*2': 16}) - self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_3, rmg_reactions)) - - - def test_map_arc_rmg_species(self): - """Test the map_arc_rmg_species() function.""" - r_map, p_map = map_arc_rmg_species(arc_reaction=ARCReaction(r_species=[ARCSpecies(label='CCjC', smiles='C[CH]C')], - p_species=[ARCSpecies(label='CjCC', smiles='[CH2]CC')]), - rmg_reaction=Reaction(reactants=[Species(smiles='C[CH]C')], - products=[Species(smiles='[CH2]CC')]), - concatenate=False) - self.assertEqual(r_map, {0: 0}) - self.assertEqual(p_map, {0: 0}) - - r_map, p_map = map_arc_rmg_species(arc_reaction=ARCReaction(r_species=[ARCSpecies(label='CCjC', smiles='C[CH]C')], - p_species=[ARCSpecies(label='CjCC', smiles='[CH2]CC')]), - rmg_reaction=Reaction(reactants=[Species(smiles='C[CH]C')], - products=[Species(smiles='[CH2]CC')])) - self.assertEqual(r_map, {0: [0]}) - self.assertEqual(p_map, {0: [0]}) - - r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_1, rmg_reaction=self.rmg_reaction_1) - self.assertEqual(r_map, {0: [0], 1: [1]}) - self.assertEqual(p_map, {0: [0], 1: [1]}) - - r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_1, rmg_reaction=self.rmg_reaction_2) - self.assertEqual(r_map, {0: [1], 1: [0]}) - self.assertEqual(p_map, {0: [0], 1: [1]}) - - r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_3, rmg_reaction=self.rmg_reaction_3) - self.assertEqual(r_map, {0: [0, 1], 1: [0, 1]}) - self.assertEqual(p_map, {0: [0]}) - - rmg_reaction_1 = Reaction(reactants=[Species(smiles='N'), Species(smiles='[H]')], - products=[Species(smiles='[NH2]'), Species(smiles='[H][H]')]) - rmg_reaction_2 = Reaction(reactants=[Species(smiles='[H]'), Species(smiles='N')], - products=[Species(smiles='[H][H]'), Species(smiles='[NH2]')]) - rmg_reaction_3 = Reaction(reactants=[Species(smiles='N'), Species(smiles='[H]')], - products=[Species(smiles='[H][H]'), Species(smiles='[NH2]')]) - arc_reaction = ARCReaction(r_species=[ARCSpecies(label='NH3', smiles='N'), ARCSpecies(label='H', smiles='[H]')], - p_species=[ARCSpecies(label='NH2', smiles='[NH2]'), ARCSpecies(label='H2', smiles='[H][H]')]) - - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_1, arc_reaction=arc_reaction) - self.assertEqual(r_map, {0: [0], 1: [1]}) - self.assertEqual(p_map, {0: [0], 1: [1]}) - - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_2, arc_reaction=arc_reaction) - self.assertEqual(r_map, {0: [1], 1: [0]}) - self.assertEqual(p_map, {0: [1], 1: [0]}) - - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_3, arc_reaction=arc_reaction) - self.assertEqual(r_map, {0: [0], 1: [1]}) - self.assertEqual(p_map, {0: [1], 1: [0]}) - rmg_reaction = Reaction(reactants=[Species(smiles='[CH3]'), Species(smiles='[CH3]')], - products=[Species(smiles='CC')]) - arc_reaction = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), - ARCSpecies(label='CH3', smiles='[CH3]')], - p_species=[ARCSpecies(label='C2H6', smiles='CC')]) - - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, arc_reaction=arc_reaction) - self.assertEqual(r_map, {0: [0, 1], 1: [0, 1]}) - self.assertEqual(p_map, {0: [0]}) - - rmg_reaction = Reaction(reactants=[Species(smiles='c1ccccc1CCC=C'), Species(smiles='CCO[O]')], - products=[Species(smiles='c1ccccc1[CH]CC=C'), Species(smiles='CCOO')]) - arc_reaction = ARCReaction(r_species=[ARCSpecies(label='butenylbenzene', smiles='c1ccccc1CCC=C'), - ARCSpecies(label='peroxyl', smiles='CCO[O]')], - p_species=[ARCSpecies(label='peroxide', smiles='CCOO'), - ARCSpecies(label='butenylbenzene_rad', smiles='c1ccccc1[CH]CC=C')]) - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, - arc_reaction=arc_reaction, - concatenate=False, - ) - self.assertEqual(r_map, {0: 0, 1: 1}) - self.assertEqual(p_map, {0: 1, 1: 0}) - - rmg_reaction = Reaction(reactants=[Species(smiles='[N]N'), Species(smiles='NN')], - products=[Species(smiles='[NH]N'), Species(smiles='[NH]N')]) - arc_reaction = ARCReaction(r_species=[ARCSpecies(label='H2NN(T)', smiles='[N]N'), - ARCSpecies(label='N2H4', smiles='NN')], - p_species=[ARCSpecies(label='N2H3', smiles='[NH]N'), - ARCSpecies(label='N2H3', smiles='[NH]N')]) - r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, - arc_reaction=arc_reaction, - concatenate=False, - ) - self.assertEqual(r_map, {0: 0, 1: 1}) - self.assertEqual(p_map, {0: 0, 1: 1}) - + # def test_get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(self): + # """Test the get_atom_indices_of_labeled_atoms_in_an_rmg_reaction() function.""" + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_1) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_1, + # rmg_reaction=rmg_reactions[0]) + # self.assertEqual(r_dict['*1'], 0) + # self.assertIn(r_dict['*2'], [1, 2, 3, 4]) + # self.assertEqual(r_dict['*3'], 5) + # self.assertEqual(p_dict['*1'], 0) + # self.assertIn(p_dict['*2'], [5, 6]) + # self.assertEqual(p_dict['*3'], 4) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_1, rmg_reactions)) + # + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_2) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_2, + # rmg_reaction=rmg_reactions[0]) + # self.assertIn(r_dict['*1'], [0, 2]) + # self.assertIn(r_dict['*2'], [3, 4, 5, 8, 9, 10]) + # self.assertEqual(r_dict['*3'], 11) + # self.assertEqual(p_dict['*1'], 0) + # self.assertIn(p_dict['*2'], [11, 12, 13]) + # self.assertEqual(p_dict['*3'], 10) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_2, rmg_reactions)) + # + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.arc_reaction_4) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.arc_reaction_4, + # rmg_reaction=rmg_reactions[0]) + # self.assertEqual(r_dict['*1'], 0) + # self.assertEqual(r_dict['*2'], 2) + # self.assertIn(r_dict['*3'], [7, 8]) + # self.assertEqual(p_dict['*1'], 0) + # self.assertEqual(p_dict['*2'], 2) + # self.assertIn(p_dict['*3'], [3, 4, 5]) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.arc_reaction_4, rmg_reactions)) + # + # for atom, symbol in zip(self.rxn_2a.r_species[0].mol.atoms, ['C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H']): + # self.assertEqual(atom.symbol, symbol) + # self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[0].radical_electrons, 0) + # self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[1].radical_electrons, 1) + # self.assertEqual(self.rxn_2a.r_species[0].mol.atoms[2].radical_electrons, 0) + # self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[0].radical_electrons, 0) + # self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[1].radical_electrons, 0) + # self.assertEqual(self.rxn_2a.p_species[0].mol.atoms[2].radical_electrons, 1) + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.rxn_2a) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.rxn_2a, + # rmg_reaction=rmg_reactions[0]) + # self.assertEqual(r_dict['*1'], 1) + # self.assertIn(r_dict['*2'], [0, 2]) + # self.assertIn(r_dict['*3'], [4, 5, 6, 7, 8, 9]) + # self.assertEqual(p_dict['*1'], 1) + # self.assertEqual(p_dict['*2'], 2) + # self.assertIn(p_dict['*3'], [3, 6]) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.rxn_2a, rmg_reactions)) + # + # for atom, symbol in zip(self.rxn_2b.r_species[0].mol.atoms, ['C', 'C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H']): + # self.assertEqual(atom.symbol, symbol) + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=self.rxn_2b) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=self.rxn_2b, + # rmg_reaction=rmg_reactions[0]) + # self.assertEqual(r_dict['*1'], 1) + # self.assertIn(r_dict['*2'], [0, 6]) + # self.assertIn(r_dict['*3'], [3, 4, 5, 7, 8, 9]) + # self.assertEqual(p_dict['*1'], 1) + # self.assertEqual(p_dict['*2'], 2) + # self.assertIn(p_dict['*3'], [3, 6]) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(self.rxn_2b, rmg_reactions)) + # + # # C3H6O + C4H9O <=> C3H5O + C4H10O + # r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O') + # r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO') + # p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O') + # p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO') + # rxn_1 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], + # r_species=[r_1, r_2], p_species=[p_1, p_2]) + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_1) + # for rmg_reaction in rmg_reactions: + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_1, + # rmg_reaction=rmg_reaction) + # for d in [r_dict, p_dict]: + # self.assertEqual(len(list(d.keys())), 3) + # keys = list(d.keys()) + # for label in ['*1', '*2', '*3']: + # self.assertIn(label, keys) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_1, rmg_reactions)) + # + # p_1 = ARCSpecies(label='C3H5O', smiles='CC=C[O]') # Use a wrong resonance structure and repeat the above. + # rxn_2 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], + # r_species=[r_1, r_2], p_species=[p_1, p_2]) + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_2) + # for rmg_reaction in rmg_reactions: + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_2, + # rmg_reaction=rmg_reaction) + # for d in [r_dict, p_dict]: + # self.assertEqual(len(list(d.keys())), 3) + # keys = list(d.keys()) + # for label in ['*1', '*2', '*3']: + # self.assertIn(label, keys) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_2, rmg_reactions)) + # + # # C3H6O + C4H9O <=> C3H5O + C4H10O + # r_1 = ARCSpecies(label='C3H6O', smiles='CCC=O', xyz=self.c3h6o_xyz) + # r_2 = ARCSpecies(label='C4H9O', smiles='[CH2]C(C)CO', xyz=self.c4h9o_xyz) + # p_1 = ARCSpecies(label='C3H5O', smiles='C[CH]C=O', xyz=self.c3h5o_xyz) + # p_2 = ARCSpecies(label='C4H10O', smiles='CC(C)CO', xyz=self.c4h10o_xyz) + # rxn_3 = ARCReaction(reactants=['C3H6O', 'C4H9O'], products=['C3H5O', 'C4H10O'], + # r_species=[r_1, r_2], p_species=[p_1, p_2]) + # rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_3) + # r_dict, p_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_3, + # rmg_reaction=rmg_reactions[0]) + # self.assertEqual(r_dict, {'*3': 10, '*1': 1, '*2': 7}) + # self.assertEqual(p_dict, {'*1': 1, '*3': 9, '*2': 16}) + # self.assertTrue(_check_r_n_p_symbols_between_rmg_and_arc_rxns(rxn_3, rmg_reactions)) + # + # def test_map_arc_rmg_species(self): + # """Test the map_arc_rmg_species() function.""" + # r_map, p_map = map_arc_rmg_species( + # arc_reaction=ARCReaction(r_species=[ARCSpecies(label='CCjC', smiles='C[CH]C')], + # p_species=[ARCSpecies(label='CjCC', smiles='[CH2]CC')]), + # rmg_reaction=Reaction(reactants=[Species(smiles='C[CH]C')], + # products=[Species(smiles='[CH2]CC')]), + # concatenate=False) + # self.assertEqual(r_map, {0: 0}) + # self.assertEqual(p_map, {0: 0}) + # + # r_map, p_map = map_arc_rmg_species( + # arc_reaction=ARCReaction(r_species=[ARCSpecies(label='CCjC', smiles='C[CH]C')], + # p_species=[ARCSpecies(label='CjCC', smiles='[CH2]CC')]), + # rmg_reaction=Reaction(reactants=[Species(smiles='C[CH]C')], + # products=[Species(smiles='[CH2]CC')])) + # self.assertEqual(r_map, {0: [0]}) + # self.assertEqual(p_map, {0: [0]}) + # + # r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_1, rmg_reaction=self.rmg_reaction_1) + # self.assertEqual(r_map, {0: [0], 1: [1]}) + # self.assertEqual(p_map, {0: [0], 1: [1]}) + # + # r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_1, rmg_reaction=self.rmg_reaction_2) + # self.assertEqual(r_map, {0: [1], 1: [0]}) + # self.assertEqual(p_map, {0: [0], 1: [1]}) + # + # r_map, p_map = map_arc_rmg_species(arc_reaction=self.arc_reaction_3, rmg_reaction=self.rmg_reaction_3) + # self.assertEqual(r_map, {0: [0, 1], 1: [0, 1]}) + # self.assertEqual(p_map, {0: [0]}) + # + # rmg_reaction_1 = Reaction(reactants=[Species(smiles='N'), Species(smiles='[H]')], + # products=[Species(smiles='[NH2]'), Species(smiles='[H][H]')]) + # rmg_reaction_2 = Reaction(reactants=[Species(smiles='[H]'), Species(smiles='N')], + # products=[Species(smiles='[H][H]'), Species(smiles='[NH2]')]) + # rmg_reaction_3 = Reaction(reactants=[Species(smiles='N'), Species(smiles='[H]')], + # products=[Species(smiles='[H][H]'), Species(smiles='[NH2]')]) + # arc_reaction = ARCReaction(r_species=[ARCSpecies(label='NH3', smiles='N'), ARCSpecies(label='H', smiles='[H]')], + # p_species=[ARCSpecies(label='NH2', smiles='[NH2]'), ARCSpecies(label='H2', smiles='[H][H]')]) + # + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_1, arc_reaction=arc_reaction) + # self.assertEqual(r_map, {0: [0], 1: [1]}) + # self.assertEqual(p_map, {0: [0], 1: [1]}) + # + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_2, arc_reaction=arc_reaction) + # self.assertEqual(r_map, {0: [1], 1: [0]}) + # self.assertEqual(p_map, {0: [1], 1: [0]}) + # + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction_3, arc_reaction=arc_reaction) + # self.assertEqual(r_map, {0: [0], 1: [1]}) + # self.assertEqual(p_map, {0: [1], 1: [0]}) + # rmg_reaction = Reaction(reactants=[Species(smiles='[CH3]'), Species(smiles='[CH3]')], + # products=[Species(smiles='CC')]) + # arc_reaction = ARCReaction(r_species=[ARCSpecies(label='CH3', smiles='[CH3]'), + # ARCSpecies(label='CH3', smiles='[CH3]')], + # p_species=[ARCSpecies(label='C2H6', smiles='CC')]) + # + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, arc_reaction=arc_reaction) + # self.assertEqual(r_map, {0: [0, 1], 1: [0, 1]}) + # self.assertEqual(p_map, {0: [0]}) + # + # rmg_reaction = Reaction(reactants=[Species(smiles='c1ccccc1CCC=C'), Species(smiles='CCO[O]')], + # products=[Species(smiles='c1ccccc1[CH]CC=C'), Species(smiles='CCOO')]) + # arc_reaction = ARCReaction(r_species=[ARCSpecies(label='butenylbenzene', smiles='c1ccccc1CCC=C'), + # ARCSpecies(label='peroxyl', smiles='CCO[O]')], + # p_species=[ARCSpecies(label='peroxide', smiles='CCOO'), + # ARCSpecies(label='butenylbenzene_rad', smiles='c1ccccc1[CH]CC=C')]) + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, + # arc_reaction=arc_reaction, + # concatenate=False, + # ) + # self.assertEqual(r_map, {0: 0, 1: 1}) + # self.assertEqual(p_map, {0: 1, 1: 0}) + # + # rmg_reaction = Reaction(reactants=[Species(smiles='[N]N'), Species(smiles='NN')], + # products=[Species(smiles='[NH]N'), Species(smiles='[NH]N')]) + # arc_reaction = ARCReaction(r_species=[ARCSpecies(label='H2NN(T)', smiles='[N]N'), + # ARCSpecies(label='N2H4', smiles='NN')], + # p_species=[ARCSpecies(label='N2H3', smiles='[NH]N'), + # ARCSpecies(label='N2H3', smiles='[NH]N')]) + # r_map, p_map = map_arc_rmg_species(rmg_reaction=rmg_reaction, + # arc_reaction=arc_reaction, + # concatenate=False, + # ) + # self.assertEqual(r_map, {0: 0, 1: 1}) + # self.assertEqual(p_map, {0: 0, 1: 1}) def test_find_equivalent_atoms_in_reactants_and_products(self): """Test the find_equivalent_atoms_in_reactants_and_products() function""" @@ -1141,7 +1075,6 @@ def test_check_species_before_mapping(self): verbose=True)) self.assertTrue(check_species_before_mapping(spc_1=self.spc1, spc_2=self.spc2, verbose=True)) - def test_get_bonds_dict(self): """Test the get_bonds_dict function.""" bond_dict = get_bonds_dict(spc=ARCSpecies(label='CH4', smiles='C')) @@ -1387,22 +1320,20 @@ def test_make_bond_changes(self): label_species_atoms([spc1]) label_species_atoms([spc2]) rxn = ARCReaction(r_species = [spc1], p_species=[spc2]) - rxn.determine_family(self.db) 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")]) - rxn.determine_family(self.db) - 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 = 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) 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_bdes(rxn, r_label_dict, True), find_all_bdes(rxn, p_label_dict, False) + r_bdes, _ = find_all_bdes(rxn, True), find_all_bdes(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, @@ -1412,15 +1343,15 @@ def test_make_bond_changes(self): def test_update_xyz(self): """tests the update_xyz function""" - spc = ARCSpecies(label="test_UX",smiles = "BrC(F)Cl") + spc = ARCSpecies(label="test_UX", smiles="BrC(F)Cl") shuffle(spc.mol.atoms) update_xyz([spc]) xyz = spc.get_xyz()["symbols"] atoms = [atom.element.symbol for atom in spc.mol.atoms] for label1, label2 in zip(atoms, xyz): self.assertEqual(label1, label2) - - spc = ARCSpecies(label="test_UX", smiles = "OCl") + + spc = ARCSpecies(label="test_UX", smiles="OCl") shuffle(spc.mol.atoms) update_xyz([spc]) xyz = spc.get_xyz()["symbols"] @@ -1428,36 +1359,36 @@ def test_update_xyz(self): for label1, label2 in zip(atoms, xyz): self.assertEqual(label1, label2) - spc = ARCSpecies(label="test_UX", smiles = "BrOCl") + spc = ARCSpecies(label="test_UX", smiles="BrOCl") shuffle(spc.mol.atoms) update_xyz([spc]) xyz = spc.get_xyz()["symbols"] atoms = [atom.element.symbol for atom in spc.mol.atoms] - for label1,label2 in zip(atoms, xyz): + for label1, label2 in zip(atoms, xyz): self.assertEqual(label1, label2) def test_add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(self): "test the add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion function" - spc1 = ARCSpecies(label="c4", smiles = "CCCC") + spc1 = ARCSpecies(label="c4", smiles="CCCC") spc2 = spc1.copy() out_dict = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(spc1, spc2, spc1.mol.atoms[0], spc2.mol.atoms[0], - [0,1,2,3], - {0: 0,1: 1,2: 2,3: 3}, + [0, 1, 2, 3], + {0: 0, 1: 1, 2: 2, 3: 3}, True) for key in range(7): self.assertEqual(out_dict[key], key) - - spc1 = ARCSpecies(label="c4", smiles = "cccc") + + spc1 = ARCSpecies(label="c4", smiles="cccc") spc2 = spc1.copy() out_dict = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(spc1, spc2, spc1.mol.atoms[0], spc2.mol.atoms[0], - [0,1,2,3], - {0: 0,1: 1,2: 2,3: 3}, + [0, 1, 2, 3], + {0: 0, 1: 1, 2: 2, 3: 3}, True) for key in range(4): self.assertEqual(out_dict[key], key) @@ -1468,8 +1399,8 @@ def test_add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(self): spc2, spc1.mol.atoms[0], spc2.mol.atoms[0], - [0,1,2,3], - {0: 0,1: 1,2: 2,3: 3}, + [0, 1, 2, 3], + {0: 0, 1: 1, 2: 2, 3: 3}, True) for key in range(4): self.assertEqual(out_dict[key], key) @@ -1477,9 +1408,8 @@ def test_add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(self): def test_find_all_bdes(self): """tests the find_all_bdes function""" rxn = ARCReaction(r_species=[self.r_1, self.r_2], p_species=[self.p_1, self.p_2]) - rxn.determine_family(self.db) - bdes = find_all_bdes(rxn, self.r_label_dict_rxn_1, True) - self.assertEqual(bdes, [(2, 4)]) + bdes = find_all_bdes(rxn, True) + self.assertEqual(bdes, [(4, 6)]) def test_copy_species_list_for_mapping(self): """tests the copy_species_list_for_mapping function""" @@ -1500,5 +1430,6 @@ def test_determine_bdes_on_spc_based_on_atom_labels(self): 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))) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) \ No newline at end of file diff --git a/arc/reaction/__init__.py b/arc/reaction/__init__.py new file mode 100644 index 0000000000..e33a984d39 --- /dev/null +++ b/arc/reaction/__init__.py @@ -0,0 +1,2 @@ +import arc.reaction.family +from arc.reaction.reaction import ARCReaction diff --git a/arc/reaction/family.py b/arc/reaction/family.py new file mode 100644 index 0000000000..090586bdb0 --- /dev/null +++ b/arc/reaction/family.py @@ -0,0 +1,827 @@ +""" +A module for working with RMG reaction families. +""" + +from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union +import ast +import os +import re + +from rmgpy.molecule import Bond, Group, Molecule + +from arc.common import clean_text, generate_resonance_structures, get_logger +from arc.imports import settings + +if TYPE_CHECKING: + from arc.species import ARCSpecies + from arc.reaction.reaction import ARCReaction + +RMG_DB_PATH = settings['RMG_DB_PATH'] +ARC_FAMILIES_PATH = settings['ARC_FAMILIES_PATH'] + + +logger = get_logger() + + +class ReactionFamily(object): + """ + A class for representing a reaction family. + + Args: + label (str): The reaction family label. + consider_arc_families (bool, optional): Whether to consider ARC's custom families + when searching for the family groups file. + + Attributes: + label (str): The reaction family label. + """ + + def __init__(self, + label: str, + consider_arc_families: bool = True, + ): + if label is None: + raise ValueError('Cannot initialize a ReactionFamily object without a label') + self.label = label + self.groups_as_lines = self.get_groups_file_as_lines(consider_arc_families=consider_arc_families) + self.reversible = is_reversible(self.groups_as_lines) + self.own_reverse = is_own_reverse(self.groups_as_lines) + self.reactants = get_reactant_groups_from_template(self.groups_as_lines) + self.reactant_num = self.get_reactant_num() + self.product_num = get_product_num(self.groups_as_lines) + entry_labels = list() + for reactant_group in self.reactants: + entry_labels.extend(reactant_group) + self.entries = get_entries(self.groups_as_lines, entry_labels=entry_labels) + self.actions = get_recipe_actions(self.groups_as_lines) + + def __str__(self): + """ + A string representation of the object. + """ + return f'ReactionFamily(label={self.label})' + + def get_groups_file_as_lines(self, consider_arc_families: bool = True) -> List[str]: + """ + Get the groups file as a list of lines. + Precedence is given to RMG families (ARC families should therefore have distinct names than RMG's) + + Args: + consider_arc_families (bool, optional): Whether to consider ARC's custom families. + + Returns: + List[str]: The groups file as a list of lines. + """ + groups_path = os.path.join(RMG_DB_PATH, 'input', 'kinetics', 'families', self.label, 'groups.py') + print(f"Looking for groups.py in RMG database at: {groups_path}") # Debug statement + if consider_arc_families: + groups_path = os.path.join(ARC_FAMILIES_PATH, f'{self.label}.py') + if not os.path.isfile(groups_path): + raise FileNotFoundError(f'Could not find the groups file for family {self.label}') + with open(groups_path, 'r') as f: + groups_as_lines = f.readlines() + return groups_as_lines + + def generate_products(self, + reactants: List['ARCSpecies'], + ) -> Dict[Union[str, Tuple[str, str]], List[Tuple[List['Molecule'], Dict[int, str]]]]: + """ + Generate a list of all the possible reaction products of this family starting from the list of ``reactants``. + + reactant_to_group_maps has the following structure:: + + {0: [{'group': 1, 'subgroup':