From 297ec6a07962c8c859c021a20236758837caeeca Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 3 Aug 2023 10:53:59 +0300 Subject: [PATCH 01/13] Minor: Simplified Species.populate_ts_checks() --- arc/species/species.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/arc/species/species.py b/arc/species/species.py index 66927c5337..99b1213bd2 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -1985,13 +1985,9 @@ def _scissors(self, def populate_ts_checks(self): """Populate (or restart) the .ts_checks attribute with default (``None``) values.""" if self.is_ts: - self.ts_checks = {'E0': None, - 'e_elect': None, - 'IRC': None, - 'freq': None, - 'normal_mode_displacement': None, - 'warnings': '', - } + keys = ['E0', 'e_elect', 'IRC', 'freq', 'normal_mode_displacement'] + self.ts_checks = {key: None for key in keys} + self.ts_checks['warnings'] = '' class TSGuess(object): From 2f6a6f87e95907e3ac14d28988a15fbe38daf341 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 24 Aug 2023 03:37:03 +0300 Subject: [PATCH 02/13] Refactor check_rxn_e_elect() and consolidate rxn energy report --- arc/checks/ts.py | 115 +++++++++++++++++++++++++---------------------- 1 file changed, 62 insertions(+), 53 deletions(-) diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 04bd3fc7be..72f0feb95b 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -64,7 +64,7 @@ def check_ts(reaction: 'ARCReaction', raise ValueError(f"Requested checks could be 'energy', 'freq', 'IRC', or 'rotors', got:\n{checks}") if 'energy' in checks or not reaction.ts_species.ts_checks['e_elect']: - check_ts_energy(reaction=reaction, verbose=verbose) + check_rxn_e_elect(reaction=reaction, verbose=verbose) if 'freq' in checks or (not reaction.ts_species.ts_checks['normal_mode_displacement'] and job is not None): check_normal_mode_displacement(reaction, job=job) @@ -99,9 +99,9 @@ def ts_passed_all_checks(species: 'ARCSpecies', return True -def check_ts_energy(reaction: 'ARCReaction', - verbose: bool = True, - ) -> None: +def check_rxn_e_elect(reaction: 'ARCReaction', + verbose: bool = True, + ) -> None: """ Check that the TS electronic energy is above both reactant and product wells in a ``reaction``. Sets the respective energy parameter 'e_elect' in the ``TS.ts_checks`` dictionary. @@ -110,45 +110,29 @@ def check_ts_energy(reaction: 'ARCReaction', reaction (ARCReaction): The reaction for which the TS is checked. verbose (bool, optional): Whether to print logging messages. """ - # Check whether E0 values are already known, e.g. from Arkane species YAML files - check_rxn_e0(reaction=reaction) + check_rxn_e0(reaction=reaction, verbose=verbose) if reaction.ts_species.ts_checks['E0']: return - - r_e_elect = None if any([spc.e_elect is None for spc in reaction.r_species]) \ - else sum(spc.e_elect * reaction.get_species_count(species=spc, well=0) for spc in reaction.r_species) - p_e_elect = None if any([spc.e_elect is None for spc in reaction.p_species]) \ - else sum(spc.e_elect * reaction.get_species_count(species=spc, well=1) for spc in reaction.p_species) - ts_e_elect = reaction.ts_species.e_elect - - if verbose and all([val is not None for val in [r_e_elect, p_e_elect, ts_e_elect]]): - min_e = extremum_list([r_e_elect, p_e_elect, ts_e_elect], return_min=True) - r_text = f'{r_e_elect - min_e:.2f} kJ/mol' if r_e_elect is not None else 'None' - ts_text = f'{ts_e_elect - min_e:.2f} kJ/mol' if ts_e_elect is not None else 'None' - p_text = f'{p_e_elect - min_e:.2f} kJ/mol' if p_e_elect is not None else 'None' - logger.info( - f'\nReaction {reaction.label} (TS {reaction.ts_label}, TSG {reaction.ts_species.chosen_ts}) ' - f'has the following path electronic energy:\n' - f'Reactants: {r_text}\n' - f'TS: {ts_text}\n' - f'Products: {p_text}') - - if all([val is not None for val in [r_e_elect, p_e_elect, ts_e_elect]]): - # We have all params, we can make a quantitative decision. - if ts_e_elect > r_e_elect + 1.0 and ts_e_elect > p_e_elect + 1.0: - # TS is above both wells. - reaction.ts_species.ts_checks['e_elect'] = True - return - # TS is not above both wells. - if verbose: - logger.error(f'TS of reaction {reaction.label} has a lower electronic energy value than expected.') - reaction.ts_species.ts_checks['e_elect'] = False - return - # We don't have any params (some are ``None``) + r_ee = sum_list_entries([r.e_elect for r in reaction.r_species], + multipliers=[reaction.get_species_count(species=r, well=0) for r in reaction.r_species]) + p_ee = sum_list_entries([p.e_elect for p in reaction.p_species], + multipliers=[reaction.get_species_count(species=p, well=1) for p in reaction.p_species]) + ts_ee = reaction.ts_species.e_elect + if verbose: + report_ts_and_wells_energy(r_e=r_ee, p_e=p_ee, ts_e=ts_ee, rxn_label=reaction.label, + ts_label=reaction.ts_label, chosen_ts=reaction.ts_species.chosen_ts, + energy_type='electronic energy') + if all([val is not None for val in [r_ee, p_ee, ts_ee]]): + if ts_ee > r_ee + 1.0 and ts_ee > p_ee + 1.0: + reaction.ts_species.ts_checks['e_elect'] = True + return + if verbose: + logger.error(f'TS of reaction {reaction.label} has a lower electronic energy value than expected.') + reaction.ts_species.ts_checks['e_elect'] = False + return if verbose: logger.info('\n') logger.warning(f"Could not get electronic energy for all species in reaction {reaction.label}.\n") - # We don't really know. reaction.ts_species.ts_checks['e_elect'] = None if 'Could not determine TS e_elect relative to the wells; ' not in reaction.ts_species.ts_checks['warnings']: reaction.ts_species.ts_checks['warnings'] += 'Could not determine TS e_elect relative to the wells; ' @@ -161,6 +145,7 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', output: dict, sp_level: 'Level', freq_scale_factor: float = 1.0, + verbose: bool = True, ) -> Optional[bool]: """ Checking the E0 values between wells and a TS in a ``reaction`` using ZPE from statmech. @@ -173,6 +158,7 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', output (dict): The Scheduler output dictionary. sp_level (Level): The single-point energy level of theory. freq_scale_factor (float, optional): The frequency scaling factor. + verbose (bool, optional): Whether to print logging messages. Returns: Optional[bool]: Whether the test failed and the scheduler should switch to a different TS guess, @@ -201,12 +187,12 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', e0_only=True, skip_rotors=True, ) - e0_pass = check_rxn_e0(reaction=rxn_copy) + e0_pass = check_rxn_e0(reaction=rxn_copy, verbose=verbose) if not e0_pass: if rxn_copy.ts_species.ts_guesses_exhausted: return False - return True # Switch TS. - return False # Don't switch TS. + return True + return False def check_rxn_e0(reaction: 'ARCReaction', @@ -229,18 +215,9 @@ def check_rxn_e0(reaction: 'ARCReaction', p_e0 = sum_list_entries([p.e0 for p in reaction.p_species], multipliers=[reaction.get_species_count(species=p, well=1) for p in reaction.p_species]) ts_e0 = reaction.ts_species.e0 - - if verbose and all([val is not None for val in [r_e0, p_e0, ts_e0]]): - min_e0 = extremum_list([r_e0, p_e0, ts_e0], return_min=True) - r_text = f'{r_e0 - min_e0:.2f} kJ/mol' if r_e0 is not None else 'None' - ts_text = f'{ts_e0 - min_e0:.2f} kJ/mol' if ts_e0 is not None else 'None' - p_text = f'{p_e0 - min_e0:.2f} kJ/mol' if p_e0 is not None else 'None' - logger.info( - f'\nReaction {reaction.label} (TS {reaction.ts_label}, TSG {reaction.ts_species.chosen_ts}) ' - f'has the following path E0 values:\n' - f'Reactants: {r_text}\n' - f'TS: {ts_text}\n' - f'Products: {p_text}') + if verbose: + report_ts_and_wells_energy(r_e=r_e0, p_e=p_e0, ts_e=ts_e0, rxn_label=reaction.label, ts_label=reaction.ts_label, + chosen_ts=reaction.ts_species.chosen_ts, energy_type='E0') if any(e0 is None for e0 in [r_e0, p_e0, ts_e0]): return None if r_e0 >= ts_e0 or p_e0 >= ts_e0: @@ -250,6 +227,38 @@ def check_rxn_e0(reaction: 'ARCReaction', return True +def report_ts_and_wells_energy(r_e: float, + p_e: float, + ts_e: float, + rxn_label: str, + ts_label: str, + chosen_ts: int, + energy_type: str = 'electronic energy', + ): + """ + Report the relative R/TS/P energies. + + Args: + r_e (float): The reactant energy. + p_e (float): The product energy. + ts_e (float): The TS energy. + rxn_label (str): The reaction label. + ts_label (str): The TS label. + chosen_ts (int): The TSG number. + energy_type (str): The energy type: 'electronic energy' or 'E0'. + """ + if all([val is not None for val in [r_e, p_e, ts_e]]): + min_e = extremum_list(lst=[r_e, p_e, ts_e], return_min=True) + r_text = f'{r_e - min_e:.2f} kJ/mol' + ts_text = f'{ts_e - min_e:.2f} kJ/mol' + p_text = f'{p_e - min_e:.2f} kJ/mol' + logger.info( + f'\nReaction {rxn_label} (TS {ts_label}, TSG {chosen_ts}) has the following {energy_type} values:\n' + f'Reactants: {r_text}\n' + f'TS: {ts_text}\n' + f'Products: {p_text}') + + def check_normal_mode_displacement(reaction: 'ARCReaction', job: Optional['JobAdapter'], amplitudes: Optional[Union[float, List[float]]] = None, From 4f6eb21fb2d08932b3f0ec7f5022d58284a27df8 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 3 Aug 2023 09:45:20 +0300 Subject: [PATCH 03/13] Tests: Refactor check_rxn_e_elect() --- arc/checks/ts_test.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/arc/checks/ts_test.py b/arc/checks/ts_test.py index 0095f32dc2..57fe84e05a 100644 --- a/arc/checks/ts_test.py +++ b/arc/checks/ts_test.py @@ -254,41 +254,41 @@ def test_did_ts_pass_all_checks(self): self.assertTrue(ts.ts_passed_all_checks(spc, exemptions=['normal_mode_displacement', 'warnings'])) spc.ts_checks['e_elect'] = False - def test_check_ts_energy(self): - """Test the check_ts_energy() function.""" + def test_check_rxn_e_elect(self): + """Test the check_rxn_e_elect() function.""" rxn1 = ARCReaction(r_species=[ARCSpecies(label='s1', smiles='C')], p_species=[ARCSpecies(label='s2', smiles='C')]) rxn1.ts_species = ARCSpecies(label='TS', is_ts=True) # no data rxn1.ts_species.populate_ts_checks() - ts.check_ts_energy(reaction=rxn1) + ts.check_rxn_e_elect(reaction=rxn1) self.assertIsNone(rxn1.ts_species.ts_checks['e_elect']) # only E0 (incorrect) rxn1.r_species[0].e0 = 2 rxn1.p_species[0].e0 = 50 rxn1.ts_species.e0 = -100 rxn1.ts_species.populate_ts_checks() - ts.check_ts_energy(reaction=rxn1) + ts.check_rxn_e_elect(reaction=rxn1) self.assertFalse(rxn1.ts_species.ts_checks['E0']) # only E0 (partial data) rxn1.r_species[0].e0 = 2 rxn1.p_species[0].e0 = None rxn1.ts_species.e0 = -100 rxn1.ts_species.populate_ts_checks() - ts.check_ts_energy(reaction=rxn1) + ts.check_rxn_e_elect(reaction=rxn1) self.assertIsNone(rxn1.ts_species.ts_checks['e_elect']) # check e_elect (correct) rxn1.r_species[0].e_elect = 2 rxn1.p_species[0].e_elect = 50 rxn1.ts_species.e_elect = 100 rxn1.ts_species.populate_ts_checks() - ts.check_ts_energy(reaction=rxn1) + ts.check_rxn_e_elect(reaction=rxn1) self.assertTrue(rxn1.ts_species.ts_checks['e_elect']) # incorrect e_elect rxn1.r_species[0].e_elect = 2 rxn1.p_species[0].e_elect = 50 rxn1.ts_species.e_elect = -100 rxn1.ts_species.populate_ts_checks() - ts.check_ts_energy(reaction=rxn1) + ts.check_rxn_e_elect(reaction=rxn1) self.assertFalse(rxn1.ts_species.ts_checks['e_elect']) def test_determine_changing_bond(self): From 802733ebb54fb5130f3f1b7e561ade410f833323 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 22 Aug 2023 11:24:30 +0300 Subject: [PATCH 04/13] Refactor compute_rxn_e0() Now it does not run the check, and returns the reaction copy with populated E0 values --- arc/checks/ts.py | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 72f0feb95b..65e5be5c98 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -138,17 +138,17 @@ def check_rxn_e_elect(reaction: 'ARCReaction', reaction.ts_species.ts_checks['warnings'] += 'Could not determine TS e_elect relative to the wells; ' -def compute_and_check_rxn_e0(reaction: 'ARCReaction', - species_dict: dict, - project_directory: str, - kinetics_adapter: str, - output: dict, - sp_level: 'Level', - freq_scale_factor: float = 1.0, - verbose: bool = True, - ) -> Optional[bool]: +def compute_rxn_e0(reaction: 'ARCReaction', + species_dict: dict, + project_directory: str, + kinetics_adapter: str, + output: dict, + sp_level: 'Level', + freq_scale_factor: float = 1.0, + ) -> Optional['ARCReaction']: """ Checking the E0 values between wells and a TS in a ``reaction`` using ZPE from statmech. + This function computed E0 values and populates them in a copy of the given reaction instance. Args: reaction (ARCReaction): The reaction for which the TS is checked. @@ -158,11 +158,9 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', output (dict): The Scheduler output dictionary. sp_level (Level): The single-point energy level of theory. freq_scale_factor (float, optional): The frequency scaling factor. - verbose (bool, optional): Whether to print logging messages. Returns: - Optional[bool]: Whether the test failed and the scheduler should switch to a different TS guess, - ``None`` if the test couldn't be performed. + Optional['ARCReaction']: A copy of the reaction object with E0 values populated. """ for spc in reaction.r_species + reaction.p_species + [reaction.ts_species]: folder = 'rxns' if species_dict[spc.label].is_ts else 'Species' @@ -187,12 +185,7 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', e0_only=True, skip_rotors=True, ) - e0_pass = check_rxn_e0(reaction=rxn_copy, verbose=verbose) - if not e0_pass: - if rxn_copy.ts_species.ts_guesses_exhausted: - return False - return True - return False + return rxn_copy def check_rxn_e0(reaction: 'ARCReaction', From 2b852e304b120e291f0e0aa2fa5d3f2aad72c3e4 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 3 Aug 2023 10:36:07 +0300 Subject: [PATCH 05/13] Tests: Modify the compute_rxn_e0() test --- arc/checks/ts_test.py | 81 ++++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 40 deletions(-) diff --git a/arc/checks/ts_test.py b/arc/checks/ts_test.py index 57fe84e05a..2688f487f2 100644 --- a/arc/checks/ts_test.py +++ b/arc/checks/ts_test.py @@ -215,6 +215,24 @@ def setUpClass(cls): (-1.1189110146736851, -0.8090395034382198, 0.6657966599789673), (-1.1265684046717404, -0.2344009055503307, -1.0127644068816903))} + cls.species_dict_8 = {spc.label: spc for spc in cls.rxn_8.r_species + cls.rxn_8.p_species + [cls.rxn_8.ts_species]} + cls.project_directory_8 = os.path.join(ts.ARC_PATH, 'Projects', 'arc_project_for_testing_delete_after_usage5') + cls.output_dict_8 = {'iC3H7': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'iC3H7.out'), + 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'iC3H7.out'), + 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'iC3H7.out'), + 'composite': ''}, + 'convergence': True}, + 'nC3H7': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'nC3H7.out'), + 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'nC3H7.out'), + 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'nC3H7.out'), + 'composite': ''}, + 'convergence': True}, + 'TS8': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'TS_nC3H7-iC3H7.out'), + 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'TS_nC3H7-iC3H7.out'), + 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'TS_nC3H7-iC3H7.out'), + 'composite': ''}, + 'convergence': True}} + def test_check_ts(self): """Test the check_ts() function.""" self.job1.local_path_to_output_file = os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'TS_C3_intraH_8.out') @@ -309,48 +327,31 @@ def test_determine_changing_bond(self): change = ts.determine_changing_bond(bond=(0, 1), dmat_bonds_1=dmat_bonds_2, dmat_bonds_2=dmat_bonds_1) self.assertEqual(change, 'breaking') - def test_compute_and_check_rxn_e0(self): - """Test the compute_and_check_rxn_e0() function.""" - species_dict = {spc.label: spc for spc in self.rxn_8.r_species + self.rxn_8.p_species + [self.rxn_8.ts_species]} - project_directory = os.path.join(ts.ARC_PATH, 'Projects', 'arc_project_for_testing_delete_after_usage5') - output_dict = {'iC3H7': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'iC3H7.out'), - 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'iC3H7.out'), - 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'iC3H7.out'), - 'composite': '', - }, - 'convergence': True, - }, - 'nC3H7': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'nC3H7.out'), - 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'nC3H7.out'), - 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'nC3H7.out'), - 'composite': '', - }, - 'convergence': True, - }, - 'TS8': {'paths': {'freq': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'freq', 'TS_nC3H7-iC3H7.out'), - 'sp': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'TS_nC3H7-iC3H7.out'), - 'opt': os.path.join(ts.ARC_PATH, 'arc', 'testing', 'opt', 'TS_nC3H7-iC3H7.out'), - 'composite': '', - }, - 'convergence': True, - }, - } + def test_compute_rxn_e0(self): + """Test the compute_rxn_e0() function.""" for spc_label in self.rxn_8.reactants + self.rxn_8.products + [self.rxn_8.ts_label]: - folder = 'rxns' if species_dict[spc_label].is_ts else 'Species' - base_path = os.path.join(project_directory, 'output', folder, spc_label, 'geometry') + folder = 'rxns' if self.species_dict_8[spc_label].is_ts else 'Species' + base_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry') os.makedirs(base_path) - freq_path = os.path.join(project_directory, 'output', folder, spc_label, 'geometry', 'freq.out') - shutil.copy(src=output_dict[spc_label]['paths']['freq'], dst=freq_path) - - check_e0 = ts.compute_and_check_rxn_e0(reaction=self.rxn_8, - species_dict=species_dict, - project_directory=project_directory, - kinetics_adapter='arkane', - output=output_dict, - sp_level=Level(repr='uhf/3-21g'), - freq_scale_factor=1.0, - ) - self.assertFalse(check_e0) # Tests whether a TS switch is needed. The E0 check passes, thus we expect ``False``. + freq_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry', 'freq.out') + shutil.copy(src=self.output_dict_8[spc_label]['paths']['freq'], dst=freq_path) + + self.assertIsNone(self.rxn_8.r_species[0].e0) + self.assertIsNone(self.rxn_8.p_species[0].e0) + self.assertIsNone(self.rxn_8.ts_species.e0) + + rxn_copy = ts.compute_rxn_e0(reaction=self.rxn_8, + species_dict=self.species_dict_8, + project_directory=self.project_directory_8, + kinetics_adapter='arkane', + output=self.output_dict_8, + sp_level=Level(repr='uhf/3-21g'), + freq_scale_factor=1.0, + ) + + self.assertAlmostEquals(rxn_copy.r_species[0].e0, -306893.334, places=1) + self.assertAlmostEquals(rxn_copy.p_species[0].e0, -306902.397, places=1) + self.assertAlmostEquals(rxn_copy.ts_species.e0, -306655.822, places=1) def test_check_normal_mode_displacement(self): """Test the check_normal_mode_displacement() function.""" From 1f700ac2fc5448c27d03d863134888d3c877917d Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 3 Aug 2023 10:40:29 +0300 Subject: [PATCH 06/13] Modified check_rxn_e0() to change test status w/o returning a value --- arc/checks/ts.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 65e5be5c98..5de9e48c02 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -190,34 +190,31 @@ def compute_rxn_e0(reaction: 'ARCReaction', def check_rxn_e0(reaction: 'ARCReaction', verbose: bool = True, - ) -> Optional[bool]: + ): """ Check the E0 values between wells and a TS in a ``reaction``, assuming that E0 values are available. Args: reaction (ARCReaction): The reaction to consider. verbose (bool, optional): Whether to print logging messages. - - Returns: - Optional[bool]: Whether the Ts E0 is above both R and P E0 wells. """ if reaction.ts_species.ts_checks['E0']: - return True + return r_e0 = sum_list_entries([r.e0 for r in reaction.r_species], multipliers=[reaction.get_species_count(species=r, well=0) for r in reaction.r_species]) p_e0 = sum_list_entries([p.e0 for p in reaction.p_species], multipliers=[reaction.get_species_count(species=p, well=1) for p in reaction.p_species]) ts_e0 = reaction.ts_species.e0 - if verbose: - report_ts_and_wells_energy(r_e=r_e0, p_e=p_e0, ts_e=ts_e0, rxn_label=reaction.label, ts_label=reaction.ts_label, - chosen_ts=reaction.ts_species.chosen_ts, energy_type='E0') if any(e0 is None for e0 in [r_e0, p_e0, ts_e0]): - return None - if r_e0 >= ts_e0 or p_e0 >= ts_e0: - reaction.ts_species.ts_checks['E0'] = False - return False - reaction.ts_species.ts_checks['E0'] = True - return True + reaction.ts_species.ts_checks['E0'] = None + else: + if verbose: + report_ts_and_wells_energy(r_e=r_e0, p_e=p_e0, ts_e=ts_e0, rxn_label=reaction.label, ts_label=reaction.ts_label, + chosen_ts=reaction.ts_species.chosen_ts, energy_type='E0') + if r_e0 >= ts_e0 or p_e0 >= ts_e0: + reaction.ts_species.ts_checks['E0'] = False + else: + reaction.ts_species.ts_checks['E0'] = True def report_ts_and_wells_energy(r_e: float, From c3460c47137147aa4b9749888f4b4e371b7982cb Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sun, 20 Aug 2023 03:02:12 +0300 Subject: [PATCH 07/13] Tests: Added check_rxn_e0() test --- arc/checks/ts_test.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/arc/checks/ts_test.py b/arc/checks/ts_test.py index 2688f487f2..20c548450b 100644 --- a/arc/checks/ts_test.py +++ b/arc/checks/ts_test.py @@ -332,7 +332,7 @@ def test_compute_rxn_e0(self): for spc_label in self.rxn_8.reactants + self.rxn_8.products + [self.rxn_8.ts_label]: folder = 'rxns' if self.species_dict_8[spc_label].is_ts else 'Species' base_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry') - os.makedirs(base_path) + os.makedirs(base_path, exist_ok=True) freq_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry', 'freq.out') shutil.copy(src=self.output_dict_8[spc_label]['paths']['freq'], dst=freq_path) @@ -353,6 +353,27 @@ def test_compute_rxn_e0(self): self.assertAlmostEquals(rxn_copy.p_species[0].e0, -306902.397, places=1) self.assertAlmostEquals(rxn_copy.ts_species.e0, -306655.822, places=1) + def test_check_rxn_e0(self): + """Test the check_rxn_e0() function.""" + for spc_label in self.rxn_8.reactants + self.rxn_8.products + [self.rxn_8.ts_label]: + folder = 'rxns' if self.species_dict_8[spc_label].is_ts else 'Species' + base_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry') + os.makedirs(base_path, exist_ok=True) + freq_path = os.path.join(self.project_directory_8, 'output', folder, spc_label, 'geometry', 'freq.out') + shutil.copy(src=self.output_dict_8[spc_label]['paths']['freq'], dst=freq_path) + rxn_copy = ts.compute_rxn_e0(reaction=self.rxn_8, + species_dict=self.species_dict_8, + project_directory=self.project_directory_8, + kinetics_adapter='arkane', + output=self.output_dict_8, + sp_level=Level(repr='uhf/3-21g'), + freq_scale_factor=1.0, + ) + + self.assertIsNone(rxn_copy.ts_species.ts_checks['E0']) + ts.check_rxn_e0(reaction=rxn_copy, verbose=True) + self.assertTrue(rxn_copy.ts_species.ts_checks['E0']) + def test_check_normal_mode_displacement(self): """Test the check_normal_mode_displacement() function.""" self.rxn_2a.ts_species.populate_ts_checks() From 3b21f71a4f7d0d4f66a9e40222a52467f87e90e2 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sat, 5 Aug 2023 20:00:14 +0300 Subject: [PATCH 08/13] Call all TS energy tests from the same place --- arc/checks/ts.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 5de9e48c02..6d2ca6cdb9 100644 --- a/arc/checks/ts.py +++ b/arc/checks/ts.py @@ -37,10 +37,16 @@ def check_ts(reaction: 'ARCReaction', - verbose: bool = True, job: Optional['JobAdapter'] = None, checks: Optional[List[str]] = None, rxn_zone_atom_indices: Optional[List[int]] = None, + species_dict: Optional[dict] = None, + project_directory: Optional[str] = None, + kinetics_adapter: Optional[str] = None, + output: Optional[dict] = None, + sp_level: Optional['Level'] = None, + freq_scale_factor: float = 1.0, + verbose: bool = True, ): """ Check the TS in terms of energy, normal mode displacement, and IRC. @@ -52,19 +58,36 @@ def check_ts(reaction: 'ARCReaction', Args: reaction (ARCReaction): The reaction for which the TS is checked. - verbose (bool, optional): Whether to print logging messages. job (JobAdapter, optional): The frequency job object instance. checks (List[str], optional): Specific checks to run. Optional values: 'energy', 'freq', 'IRC', 'rotors'. rxn_zone_atom_indices (List[int], optional): The 0-indices of atoms identified by the normal mode displacement as the reaction zone. Automatically determined if not given. + species_dict (dict, optional): The Scheduler species dictionary. + project_directory (str, optional): The path to ARC's project directory. + kinetics_adapter (str, optional): The statmech software to use for kinetic rate coefficient calculations. + output (dict, optional): The Scheduler output dictionary. + sp_level (Level, optional): The single-point energy level of theory. + freq_scale_factor (float, optional): The frequency scaling factor. + verbose (bool, optional): Whether to print logging messages. """ checks = checks or list() for entry in checks: if entry not in ['energy', 'freq', 'IRC', 'rotors']: raise ValueError(f"Requested checks could be 'energy', 'freq', 'IRC', or 'rotors', got:\n{checks}") - if 'energy' in checks or not reaction.ts_species.ts_checks['e_elect']: - check_rxn_e_elect(reaction=reaction, verbose=verbose) + if 'energy' in checks: + if not reaction.ts_species.ts_checks['E0']: + rxn_copy = compute_rxn_e0(reaction=reaction, + species_dict=species_dict, + project_directory=project_directory, + kinetics_adapter=kinetics_adapter, + output=output, + sp_level=sp_level, + freq_scale_factor=freq_scale_factor) + reaction.copy_e0_values(rxn_copy) + check_rxn_e0(reaction=reaction, verbose=verbose) + if reaction.ts_species.ts_checks['E0'] is None and not reaction.ts_species.ts_checks['e_elect']: + check_rxn_e_elect(reaction=reaction, verbose=verbose) if 'freq' in checks or (not reaction.ts_species.ts_checks['normal_mode_displacement'] and job is not None): check_normal_mode_displacement(reaction, job=job) @@ -162,6 +185,9 @@ def compute_rxn_e0(reaction: 'ARCReaction', Returns: Optional['ARCReaction']: A copy of the reaction object with E0 values populated. """ + if any(val is None for val in [species_dict, project_directory, kinetics_adapter, + output, sp_level, freq_scale_factor]): + return None for spc in reaction.r_species + reaction.p_species + [reaction.ts_species]: folder = 'rxns' if species_dict[spc.label].is_ts else 'Species' freq_path = os.path.join(project_directory, 'output', folder, spc.label, 'geometry', 'freq.out') From d1417c6578c1b12b6e9b76f9d7c9619b8c1056a7 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sat, 5 Aug 2023 18:39:15 +0300 Subject: [PATCH 09/13] Check reaction E0 through check_ts() --- arc/scheduler.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/arc/scheduler.py b/arc/scheduler.py index 2b1662674c..7f7c033c50 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -17,7 +17,7 @@ import arc.rmgdb as rmgdb from arc import parser, plotter from arc.checks.common import get_i_from_job_name, sum_time_delta -from arc.checks.ts import check_imaginary_frequencies, check_ts, check_irc_species_and_rxn, compute_and_check_rxn_e0 +from arc.checks.ts import check_imaginary_frequencies, check_ts, check_irc_species_and_rxn from arc.common import (extremum_list, get_angle_in_180_range, get_logger, @@ -2193,8 +2193,6 @@ def parse_composite_geo(self, else: self.output[label]['isomorphism'] += 'composite did not pass isomorphism check; ' success &= is_isomorphic - if success: - self.check_rxn_e0_by_spc(label) return success elif not self.species_dict[label].is_ts and self.trsh_ess_jobs: self.troubleshoot_negative_freq(label=label, job=job) @@ -2431,20 +2429,22 @@ def check_rxn_e0_by_spc(self, label: str): """ for rxn in self.rxn_list: labels = rxn.reactants + rxn.products + [rxn.ts_label] - if label in labels and not rxn.ts_species.ts_checks['E0'] \ + if label in labels and rxn.ts_species.ts_checks['E0'] is None \ and all([(species_has_sp(output_dict) and species_has_freq(output_dict)) or self.species_dict[spc_label].yml_path is not None for spc_label, output_dict in self.output.items() if spc_label in labels]): - switch_ts = compute_and_check_rxn_e0(reaction=rxn, - species_dict=self.species_dict, - project_directory=self.project_directory, - kinetics_adapter=self.kinetics_adapter, - output=self.output, - sp_level=self.sp_level if not self.composite_method else self.composite_method, - freq_scale_factor=self.freq_scale_factor, - ) - if switch_ts is True: - logger.info(f'TS status for reaction {rxn.label} is:\n{rxn.ts_species.ts_checks}.\n' + check_ts(reaction=rxn, + checks=['energy'], + species_dict=self.species_dict, + project_directory=self.project_directory, + kinetics_adapter=self.kinetics_adapter, + output=self.output, + sp_level=self.sp_level if not self.composite_method else self.composite_method, + freq_scale_factor=self.freq_scale_factor, + verbose=True, + ) + if rxn.ts_species.ts_checks['E0'] is False: + logger.info(f'TS {rxn.ts_species.label} of reaction {rxn.label} did not pass the E0 check.\n' f'Switching TS.\n') self.switch_ts(rxn.ts_label) From 999cd6b2771620616dda71c31de84e1d6d353ad6 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Fri, 4 Aug 2023 17:26:02 +0300 Subject: [PATCH 10/13] Check species yml_path in sp/geo/freq checks in Scheduler Assuming that if the Arkane YAML path is not None, the the species has sp/geo/freq --- arc/scheduler.py | 46 +++++++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/arc/scheduler.py b/arc/scheduler.py index 7f7c033c50..75ed799dbd 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -2332,7 +2332,8 @@ def check_freq_job(self, self.output[label]['warnings'] += wrong_freq_message if job.job_status[1]['status'] != 'done' or (not freq_ok and not self.species_dict[label].is_ts): self.troubleshoot_ess(label=label, job=job, level_of_theory=job.level) - if job.job_status[1]['status'] == 'done' and freq_ok and not switch_ts and species_has_sp(self.output[label]): + if (job.job_status[1]['status'] == 'done' and freq_ok and not switch_ts + and species_has_sp(self.output[label], self.species_dict[label].yml_path)): self.check_rxn_e0_by_spc(label) def check_negative_freq(self, @@ -2430,8 +2431,8 @@ def check_rxn_e0_by_spc(self, label: str): for rxn in self.rxn_list: labels = rxn.reactants + rxn.products + [rxn.ts_label] if label in labels and rxn.ts_species.ts_checks['E0'] is None \ - and all([(species_has_sp(output_dict) and species_has_freq(output_dict)) - or self.species_dict[spc_label].yml_path is not None + and all([(species_has_sp(output_dict, self.species_dict[spc_label].yml_path) + and species_has_freq(output_dict, self.species_dict[spc_label].yml_path)) for spc_label, output_dict in self.output.items() if spc_label in labels]): check_ts(reaction=rxn, checks=['energy'], @@ -2445,7 +2446,7 @@ def check_rxn_e0_by_spc(self, label: str): ) if rxn.ts_species.ts_checks['E0'] is False: logger.info(f'TS {rxn.ts_species.label} of reaction {rxn.label} did not pass the E0 check.\n' - f'Switching TS.\n') + f'Searching for a better TS conformer...\n') self.switch_ts(rxn.ts_label) def switch_ts(self, label: str): @@ -2549,20 +2550,8 @@ def post_sp_actions(self, self.output[label]['paths']['sp_no_sol'] = sp_path self.output[label]['paths']['sp'] = original_sp_path # restore the original path - if self.species_dict[label].is_ts: - for rxn in self.rxn_dict.values(): - if rxn.ts_label == label: - if not rxn.ts_species.ts_checks['e_elect']: - check_ts(reaction=rxn, verbose=True, checks=['energy']) - if species_has_freq(self.output[label]) and not rxn.ts_species.ts_checks['E0']: - self.check_rxn_e0_by_spc(label) - if not (rxn.ts_species.ts_checks['E0'] or rxn.ts_species.ts_checks['e_elect']) \ - and (self.output[label]['paths']['freq'] or self.species_dict[label].e0): - logger.info(f'TS {label} did not pass the energy check. ' - f'Status is:\n{self.species_dict[label].ts_checks}\n' - f'Searching for a better TS conformer...') - self.switch_ts(label=label) - break + if species_has_freq(self.output[label], self.species_dict[label].yml_path): + self.check_rxn_e0_by_spc(label) # set *at the end* to differentiate between sp jobs when using complex solvation corrections self.output[label]['job_types']['sp'] = True @@ -3627,46 +3616,61 @@ def generate_final_ts_guess_report(self): save_yaml_file(path=path, content=content) -def species_has_freq(species_output_dict: dict) -> bool: +def species_has_freq(species_output_dict: dict, + yml_path: Optional[str] = None, + ) -> bool: """ Checks whether a species has valid converged frequencies using it's output dict. Args: species_output_dict (dict): The species output dict (i.e., Scheduler.output[label]). + yml_path (str): THe species Arkane YAML file path. Returns: bool Whether a species has valid converged frequencies. """ + if yml_path is not None: + return True if species_output_dict['paths']['freq'] or species_output_dict['paths']['composite']: return True return False -def species_has_geo(species_output_dict: dict) -> bool: +def species_has_geo(species_output_dict: dict, + yml_path: Optional[str] = None, + ) -> bool: """ Checks whether a species has a valid converged geometry using it's output dict. Args: species_output_dict (dict): The species output dict (i.e., Scheduler.output[label]). + yml_path (str): THe species Arkane YAML file path. Returns: bool Whether a species has a valid converged geometry. """ + if yml_path is not None: + return True if species_output_dict['paths']['geo'] or species_output_dict['paths']['composite']: return True return False -def species_has_sp(species_output_dict: dict) -> bool: +def species_has_sp(species_output_dict: dict, + yml_path: Optional[str] = None, + ) -> bool: """ Checks whether a species has a valid converged single-point energy using it's output dict. Args: species_output_dict (dict): The species output dict (i.e., Scheduler.output[label]). + yml_path (str): THe species Arkane YAML file path. Returns: bool Whether a species has a valid converged single-point energy. """ + if yml_path is not None: + return True if species_output_dict['paths']['sp'] or species_output_dict['paths']['composite']: return True return False From 397e303e62fcaf5ee590508b686a548769f59d6c Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Fri, 4 Aug 2023 17:32:14 +0300 Subject: [PATCH 11/13] Added species_has_sp_and_freq() into Scheduler --- arc/scheduler.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/arc/scheduler.py b/arc/scheduler.py index 75ed799dbd..b52b54aef6 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -2431,8 +2431,7 @@ def check_rxn_e0_by_spc(self, label: str): for rxn in self.rxn_list: labels = rxn.reactants + rxn.products + [rxn.ts_label] if label in labels and rxn.ts_species.ts_checks['E0'] is None \ - and all([(species_has_sp(output_dict, self.species_dict[spc_label].yml_path) - and species_has_freq(output_dict, self.species_dict[spc_label].yml_path)) + and all([species_has_sp_and_freq(output_dict, self.species_dict[spc_label].yml_path) for spc_label, output_dict in self.output.items() if spc_label in labels]): check_ts(reaction=rxn, checks=['energy'], @@ -3675,3 +3674,18 @@ def species_has_sp(species_output_dict: dict, return True return False + +def species_has_sp_and_freq(species_output_dict: dict, + yml_path: Optional[str] = None, + ) -> bool: + """ + Checks whether a species has a valid converged single-point energy and valid converged frequencies. + + Args: + species_output_dict (dict): The species output dict (i.e., Scheduler.output[label]). + yml_path (str): THe species Arkane YAML file path. + + Returns: bool + Whether a species has a valid converged single-point energy and frequencies. + """ + return species_has_sp(species_output_dict, yml_path) and species_has_freq(species_output_dict, yml_path) From 659bf33506a6085d4a9eb9f1aa1c93c2e6cb7ad9 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sat, 5 Aug 2023 20:00:04 +0300 Subject: [PATCH 12/13] Added copy_e0_values() method to ARCReacrtion --- arc/reaction.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/arc/reaction.py b/arc/reaction.py index f9f42637de..09e9dfbe42 100644 --- a/arc/reaction.py +++ b/arc/reaction.py @@ -960,6 +960,21 @@ def get_bonds(self) -> Tuple[list, list]: len_atoms += spc.number_of_atoms return r_bonds, p_bonds + def copy_e0_values(self, other_rxn: Optional['ARCReaction']): + """ + Copy the E0 values from another reaction object instance for the TS + and for all species if they have corresponding labels. + + Args: + other_rxn (ARCReaction): An ARCReaction object instance from which E0 values will be copied. + """ + if other_rxn is not None: + self.ts_species.e0 = self.ts_species.e0 or other_rxn.ts_species.e0 + for spc in self.r_species + self.p_species: + for other_spc in other_rxn.r_species + other_rxn.p_species: + if spc.label == other_spc.label: + spc.e0 = spc.e0 or other_spc.e0 + def remove_dup_species(species_list: List[ARCSpecies]) -> List[ARCSpecies]: """ From e8e3fe411d1e0a077d6f9de3a5d76a0b0b791a83 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 22 Aug 2023 11:24:39 +0300 Subject: [PATCH 13/13] Tests: Scheduler species_has_x() tests for species with `yml_path` attr --- arc/scheduler_test.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/arc/scheduler_test.py b/arc/scheduler_test.py index e841348c64..83d11b94f8 100644 --- a/arc/scheduler_test.py +++ b/arc/scheduler_test.py @@ -16,7 +16,7 @@ from arc.job.factory import job_factory from arc.level import Level from arc.plotter import save_conformers_file -from arc.scheduler import Scheduler, species_has_freq, species_has_geo, species_has_sp +from arc.scheduler import Scheduler, species_has_freq, species_has_geo, species_has_sp, species_has_sp_and_freq from arc.imports import settings from arc.reaction import ARCReaction from arc.species.species import ARCSpecies @@ -667,6 +667,12 @@ def test_species_has_geo_sp_freq(self): self.assertTrue(species_has_property((species_output_dict))) species_output_dict = {'paths': {property_: True, 'composite': True}} self.assertTrue(species_has_property((species_output_dict))) + yml_path=os.path.join(ARC_PATH, 'arc', 'testing', 'yml_testing', 'N4H6.yml') + species_output_dict = {'paths': {'geo': False, 'sp': False, 'freq': False, 'composite': False}} + self.assertTrue(species_has_freq(species_output_dict=species_output_dict, yml_path=yml_path)) + self.assertTrue(species_has_geo(species_output_dict=species_output_dict, yml_path=yml_path)) + self.assertTrue(species_has_sp(species_output_dict=species_output_dict, yml_path=yml_path)) + self.assertTrue(species_has_sp_and_freq(species_output_dict=species_output_dict, yml_path=yml_path)) def test_add_label_to_unique_species_labels(self): """Test the add_label_to_unique_species_labels() method."""