diff --git a/arc/checks/ts.py b/arc/checks/ts.py index 04bd3fc7be..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_ts_energy(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) @@ -99,9 +122,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,60 +133,45 @@ 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; ' -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, - ) -> 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. @@ -175,9 +183,11 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', freq_scale_factor (float, optional): The frequency scaling factor. 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. """ + 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') @@ -201,53 +211,68 @@ def compute_and_check_rxn_e0(reaction: 'ARCReaction', e0_only=True, skip_rotors=True, ) - e0_pass = check_rxn_e0(reaction=rxn_copy) - 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 rxn_copy 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 any(e0 is None for e0 in [r_e0, p_e0, ts_e0]): + 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, + 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. - 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' + 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 {reaction.label} (TS {reaction.ts_label}, TSG {reaction.ts_species.chosen_ts}) ' - f'has the following path E0 values:\n' + 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}') - 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 def check_normal_mode_displacement(reaction: 'ARCReaction', diff --git a/arc/checks/ts_test.py b/arc/checks/ts_test.py index 0095f32dc2..20c548450b 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') @@ -254,41 +272,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): @@ -309,48 +327,52 @@ 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') - 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``. + 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) + + 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_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.""" 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]: """ diff --git a/arc/scheduler.py b/arc/scheduler.py index 2b1662674c..b52b54aef6 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) @@ -2334,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, @@ -2431,21 +2430,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'] \ - and all([(species_has_sp(output_dict) and species_has_freq(output_dict)) - or self.species_dict[spc_label].yml_path is not None + if label in labels and rxn.ts_species.ts_checks['E0'] is None \ + 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]): - 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' - f'Switching TS.\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'Searching for a better TS conformer...\n') self.switch_ts(rxn.ts_label) def switch_ts(self, label: str): @@ -2549,20 +2549,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,47 +3615,77 @@ 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 + +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) 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.""" 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):