Fixes TS E0 checks (#664)

This PR primarily improves ARC E0 checks for a reaction.
Sometimes E0 values are known (e.g., from Arkane YAML files) and we
don't have the e_elect values.
ARC should use the E0 values instead when comparing relative reactants -
TS - products well heights

Other minor additions in this PR:
- Initialize TSGuess object derived from user guesses with
`success=True` (important for facilitating user TS guesses)
- Added a shortcut argument `xyz` to `ARCReaction` object (more
intuitive in input files)
- Save the species `adjlist` if given, so that ARC/RMG won't regenerate
it for libraries, and perhaps get it wrong from xyz perception

arc/checks/ts.py

@@ -110,24 +110,28 @@ 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)
+    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
-    min_e = extremum_list([r_e_elect, p_e_elect, ts_e_elect], return_min=True)
-
-    if any([val is not None for val in [r_e_elect, p_e_elect, ts_e_elect]]):
-        if verbose:
-            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 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.
@@ -150,14 +154,14 @@ def check_ts_energy(reaction: 'ARCReaction',
         reaction.ts_species.ts_checks['warnings'] += 'Could not determine TS e_elect relative to the wells; '
 
 
-def 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_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]:
     """
     Checking the E0 values between wells and a TS in a ``reaction`` using ZPE from statmech.
 
@@ -182,7 +186,7 @@ def check_rxn_e0(reaction: 'ARCReaction',
     considered_labels = list()
     rxn_copy = reaction.copy()
     for species in rxn_copy.r_species + rxn_copy.p_species + [rxn_copy.ts_species]:
-        if species.label in considered_labels:
+        if species.label in considered_labels or species.e0:
             continue
         considered_labels.append(species.label)
         statmech_adapter = statmech_factory(statmech_adapter_label=kinetics_adapter,
@@ -197,20 +201,55 @@ def check_rxn_e0(reaction: 'ARCReaction',
                                         e0_only=True,
                                         skip_rotors=True,
                                         )
-    r_e0 = sum_list_entries([r.e0 for r in rxn_copy.r_species],
-                            multipliers=[rxn_copy.get_species_count(species=r, well=0) for r in rxn_copy.r_species])
-    p_e0 = sum_list_entries([p.e0 for p in rxn_copy.p_species],
-                            multipliers=[rxn_copy.get_species_count(species=p, well=1) for p in rxn_copy.p_species])
-    if any(e0 is None for e0 in [r_e0, p_e0, rxn_copy.ts_species.e0]) \
-            or r_e0 >= rxn_copy.ts_species.e0 or p_e0 >= rxn_copy.ts_species.e0:
-        reaction.ts_species.ts_checks['E0'] = False
+    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.
-    reaction.ts_species.ts_checks['E0'] = True
     return False  # Don't switch TS.
 
 
+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
+    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 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 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',
                                    job: Optional['JobAdapter'],
                                    amplitudes: Optional[Union[float, List[float]]] = None,

arc/checks/ts_test.py

@@ -309,8 +309,8 @@ 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_check_rxn_e0(self):
-        """Test the check_rxn_e0() function."""
+    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'),
@@ -342,14 +342,14 @@ def test_check_rxn_e0(self):
             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.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,
-                                   )
+        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``.
 
     def test_check_normal_mode_displacement(self):

arc/common.py

@@ -1614,6 +1614,12 @@ def safe_copy_file(source: str,
             shutil.copyfile(src=source, dst=destination)
         except shutil.SameFileError:
             break
+        except NotADirectoryError:
+            print(f'Expected "source" and "destination" to be directories. Trying to extract base paths.')
+            if os.path.isfile(source):
+                source = os.path.dirname(source)
+            if os.path.isfile(destination):
+                destination = os.path.dirname(destination)
         except OSError:
             time.sleep(wait)
         else:

arc/main.py

@@ -227,7 +227,7 @@ def __init__(self,
                  adaptive_levels: Optional[dict] = None,
                  allow_nonisomorphic_2d: bool = False,
                  arkane_level_of_theory: Optional[Union[dict, Level, str]] = None,
-                 bac_type: str = 'p',
+                 bac_type: Optional[str] = 'p',
                  bath_gas: Optional[str] = None,
                  calc_freq_factor: bool = True,
                  compare_to_rmg: bool = True,

arc/main_test.py

@@ -16,7 +16,6 @@
 from arc.main import ARC, StatmechEnum, process_adaptive_levels
 from arc.species.species import ARCSpecies
 
-
 servers = settings['servers']
 
 
@@ -447,32 +446,32 @@ def test_process_level_of_theory(self):
         """
         Tests the process_level_of_theory function.
         """
-        arc0 = ARC(project='test_0', level_of_theory='ccsd(t)-f12/cc-pvdz-f12//b3lyp/6-311+g(3df,2p)', bac_type=False, freq_scale_factor=1,)
+        arc0 = ARC(project='test_0', level_of_theory='ccsd(t)-f12/cc-pvdz-f12//b3lyp/6-311+g(3df,2p)',
+                   bac_type=None, freq_scale_factor=1)
         arc1 = ARC(project='test_1', level_of_theory='wb97xd/6-311+g(2d,2p)',
-                                    arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
-                                    bac_type=False,
-                                    freq_scale_factor=1,
-                                    job_types= {"freq": True,
-                                                "sp"  : True,
-                                                "opt" : False,})
+                   arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                   bac_type=None,
+                   freq_scale_factor=1,
+                   job_types={"freq": True,
+                              "sp": True,
+                              "opt": False})
         arc2 = ARC(project='test_2', sp_level='wb97xd/6-311+g(2d,2p)',
-                                     opt_level='wb97xd/6-311+g(2d,2p)',
-                                     arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
-                                     bac_type=False,
-                                     freq_scale_factor=1,
-                                     job_types= {"freq": True,
-                                                 "sp"  : False,
-                                                 "opt" : False,})
+                   opt_level='wb97xd/6-311+g(2d,2p)',
+                   arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                   bac_type=None,
+                   freq_scale_factor=1,
+                   job_types={"freq": True,
+                              "sp": False,
+                              "opt": False})
         arc3 = ARC(project='test_3', sp_level='wb97xd/6-311+g(2d,2p)',
-                                     opt_level='wb97xd/6-311+g(2d,2p)',
-                                     arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
-                                     bac_type=False,
-                                     freq_scale_factor=1,
-                                     job_types= {"opt" : False,})
+                   opt_level='wb97xd/6-311+g(2d,2p)',
+                   arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                   bac_type=None,
+                   freq_scale_factor=1,
+                   job_types={"opt": False})
 
         arc0.process_level_of_theory(), arc1.process_level_of_theory(), arc2.process_level_of_theory(), arc3.process_level_of_theory()
         for arc in [arc0, arc1, arc2, arc3]:
-            print("arc")
             self.assertIsInstance(arc.sp_level, Level)
             self.assertIsInstance(arc.opt_level, Level)
             self.assertIsInstance(arc.freq_level, Level)

arc/plotter.py

@@ -777,11 +777,11 @@ def save_thermo_lib(species_list: list,
                 try:
                     thermo_library.load_entry(index=i,
                                               label=spc.label,
-                                              molecule=spc.mol_list[0].copy(deep=True).to_adjacency_list(),
+                                              molecule=spc.adjlist or spc.mol_list[0].copy(deep=True).to_adjacency_list(),
                                               thermo=spc.thermo,
                                               shortDesc=spc.thermo.comment,
                                               longDesc=spc.long_thermo_description)
-                except (InvalidAdjacencyListError, DatabaseError, ValueError) as e:
+                except (InvalidAdjacencyListError, DatabaseError, ValueError, TypeError) as e:
                     logger.error(f'Could not save species {spc.label} in the thermo library, got:')
                     logger.info(e)
             else:
@@ -806,11 +806,11 @@ def save_transport_lib(species_list, path, name, lib_long_desc=''):
                 try:
                     transport_library.load_entry(index=i,
                                                  label=spc.label,
-                                                 molecule=spc.mol_list[0].copy(deep=True).to_adjacency_list(),
+                                                 molecule=spc.adjlist or spc.mol_list[0].copy(deep=True).to_adjacency_list(),
                                                  transport=spc.transport_data,
                                                  shortDesc=spc.thermo.comment,
                                                  longDesc=description)
-                except (InvalidAdjacencyListError, ValueError) as e:
+                except (InvalidAdjacencyListError, DatabaseError, ValueError, TypeError) as e:
                     logger.error(f'Could not save species {spc.label} in the transport library, got:')
                     logger.info(e)
                 logger.info(f'\n\nTransport properties for {spc.label}:')