Resolve H_Abstraction well ordering per product dict; skip linear-segment torsions

arc/job/adapters/ts/heuristics.py

@@ -830,26 +830,33 @@ def find_distant_neighbor(mol: 'Molecule',
 
 def are_h_abs_wells_reversed(rxn: ARCReaction,
                              product_dict: dict,
-                             ) -> tuple[bool, bool]:
+                             ) -> tuple[bool, bool, bool]:
     """
     Determine whether the reactants or the products in an H_Abstraction reaction are reversed
     relative to the RMG template: R(*1)-H(*2) + R(*3)j <=> R(*1)j + R(*3)-H(*2)
+    ``reactants_reversed`` is True when the atom labeled *2 belongs to the reaction's second
+    reactant. ``products_reversed`` is True when the product bearing the atom labeled *2 is the
+    reaction's first product. ``dict_products_reversed`` is True when the atom labeled *2 belongs
+    to the first molecule of ``product_dict['products']``. ``product_dict['r_label_map']`` holds
+    zero-based global indices into the reaction's concatenated reactants, and
+    ``product_dict['p_label_map']`` holds zero-based global indices into the concatenated
+    ``product_dict['products']``.
 
     Args:
         rxn (ARCReaction): The ARCReaction object.
         product_dict (dict): The product dictionary.
 
     Returns:
-        tuple[bool, bool]: reactants_reversed, products_reversed.
+        tuple[bool, bool, bool]: reactants_reversed, products_reversed, dict_products_reversed.
     """
     r_star_2 = product_dict['r_label_map']['*2']
     p_star_2 = product_dict['p_label_map']['*2']
-    r_species, p_species = rxn.get_reactants_and_products(return_copies=True)
-    reactants_reversed = len(r_species[0].mol.atoms) < r_star_2
-    products_reversed = len(product_dict['products'][0].atoms) >= p_star_2
+    r_species, p_species = rxn.get_reactants_and_products(return_copies=False)
+    reactants_reversed = r_star_2 >= len(r_species[0].mol.atoms)
+    dict_products_reversed = p_star_2 < len(product_dict['products'][0].atoms)
     same_order_between_rxn_prods_and_dict_prods = p_species[0].is_isomorphic(product_dict['products'][0])
-    products_reversed = products_reversed == same_order_between_rxn_prods_and_dict_prods
-    return reactants_reversed, products_reversed
+    products_reversed = dict_products_reversed == same_order_between_rxn_prods_and_dict_prods
+    return reactants_reversed, products_reversed, dict_products_reversed
 
 
 def h_abstraction(reaction: ARCReaction,
@@ -860,6 +867,10 @@ def h_abstraction(reaction: ARCReaction,
                   ) -> list[dict]:
     """
     Generate TS guesses for reactions of the RMG ``H_Abstraction`` family.
+    Only the ``H_Abstraction`` entries of ``reaction.product_dicts`` are considered, and the
+    orientation of the reactant and product wells relative to the RMG template is resolved
+    separately for each of them. An entry whose *2 label does not resolve to a hydrogen atom of the
+    respective reactant is skipped.
 
     Args:
         reaction: An ARCReaction instance.
@@ -871,20 +882,26 @@ def h_abstraction(reaction: ARCReaction,
         dihedral_increment (int, optional): The dihedral increment to use for B-H-A-C and D-B-H-C dihedral scans.
 
     Returns: list[dict]
-        Entries are Cartesian coordinates of TS guesses for all reactions. Returns an empty list if
-        ``reaction.product_dicts`` is empty (i.e., no matching H_Abstraction family products were identified).
+        Entries are Cartesian coordinates of TS guesses for all reactions.
+        An empty list is returned if none of ``reaction.product_dicts`` belongs to
+        the ``H_Abstraction`` family.
     """
     xyz_guesses = list()
     dihedral_increment = dihedral_increment or DIHEDRAL_INCREMENT
-    if not reaction.product_dicts:
+    product_dicts = [product_dict for product_dict in reaction.product_dicts
+                     if product_dict.get('family') == 'H_Abstraction']
+    if not product_dicts:
         logger.warning(f'Could not generate H_Abstraction TS guesses for reaction {reaction}: '
-                        f'no product_dicts were identified for this reaction. Other TS search methods will be attempted.')
+                       f'no H_Abstraction product dictionaries were identified for it '
+                       f'(it has {len(reaction.product_dicts)} product dictionaries). '
+                       f'Other TS search methods will be attempted.')
         return xyz_guesses
-    reactants_reversed, products_reversed = are_h_abs_wells_reversed(rxn=reaction, product_dict=reaction.product_dicts[0])
-    for product_dict in reaction.product_dicts:
+    for product_dict in product_dicts:
         # Identify R1H and R2H in the "R1H + R2 <=> R1 + R2H" or "R2 + R1H <=> R2H + R1" reaction
         # The expected RMG atom labels are: R(*1)-H(*2) + R(*3)j <=> R(*1)j + R(*3)-H(*2).
         # They appear in each product_dict under the 'r_label_map' key.
+        reactants_reversed, products_reversed, dict_prods_reversed = are_h_abs_wells_reversed(
+            rxn=reaction, product_dict=product_dict)
         reactants, products = reaction.get_reactants_and_products(return_copies=False)
         reactant = reactants[int(reactants_reversed)]  # Get R(*1)-H(*2).
         reactant_2 = reactants[int(not reactants_reversed)]  # Get R(*3)j.
@@ -894,10 +911,13 @@ def h_abstraction(reaction: ARCReaction,
             # Don't modify dihedrals for an attacking H (or other linear radical) at a linear angle, C ~ A -- H1 - H2 -- H.
             dihedral_increment = 360
         h1 = product_dict['r_label_map']['*2']
-        if reactants_reversed and h1 >= len(reactants[0].mol.atoms):
+        if reactants_reversed:
             h1 -= len(reactants[0].mol.atoms)
+        if h1 >= len(r_mol.atoms) or not r_mol.atoms[h1].is_hydrogen():
+            logger.warning(f'Not using a product dictionary of reaction {reaction} in which the atom labeled *2 '
+                           f'is not a hydrogen atom of {reactant.label}.')
+            continue
         h2 = product_dict['p_label_map']['*2']
-        dict_prods_reversed = h2 < len(product_dict['products'][0].atoms)
         dict_product = product_dict['products'][int(not dict_prods_reversed)]  # Get R(*3)-H(*2) from the product_dict.
         product_atom_map = map_two_species(spc_1=dict_product, spc_2=product)
         if h2 >= len(product_atom_map):