Added a module for working with RMG reaction families (#754)

Now we can identify for each ARC reaction all possible pathways between
the reactants and products according to reaction families supported by
RMG.
Also added dummy ARC families (w/o kinetics) to assist in identifying
reaction templates.
Tests added.

Author: calvinp0

.github/workflows/cont_int.yml

@@ -36,7 +36,9 @@ jobs:
 
     - name: Cache RMG-Py
       id: cache-rmg-py
-      uses: actions/cache@v2
+      uses: actions/cache@v3
+      env:
+        CACHE_NUMBER: 2
       with:
         path: RMG-Py
         key: ${{ runner.os }}-rmg-main
@@ -54,7 +56,9 @@ jobs:
 
     - name: Cache RMG-database
       id: cache-rmg-db
-      uses: actions/cache@v2
+      uses: actions/cache@v3
+      env:
+        CACHE_NUMBER: 2
       with:
           path: RMG-database
           key: ${{ runner.os }}-rmgdb-main
@@ -71,7 +75,7 @@ jobs:
 
     - name: Cache AutoTST
       id: cache-autotst
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
           path: AutoTST
           key: ${{ runner.os }}-autotst-main
@@ -86,26 +90,26 @@ jobs:
         ref: main
         fetch-depth: 1
 
-    - name: Cache TS-GCN
-      id: cache-tsgcn
-      uses: actions/cache@v2
-      with:
-          path: TS-GCN
-          key: ${{ runner.os }}-tsgcn-main
-          restore-keys: |
-            ${{ runner.os }}-tsgcn-
-    - name: Checkout TS-GCN
-      if: steps.cache-tsgcn.outputs.cache-hit != 'true'
-      uses: actions/checkout@v3
-      with:
-        repository: ReactionMechanismGenerator/TS-GCN
-        path: TS-GCN
-        ref: main
-        fetch-depth: 1
-    
+#    - name: Cache TS-GCN
+#      id: cache-tsgcn
+#      uses: actions/cache@v3
+#      with:
+#          path: TS-GCN
+#          key: ${{ runner.os }}-tsgcn-main
+#          restore-keys: |
+#            ${{ runner.os }}-tsgcn-
+#    - name: Checkout TS-GCN
+#      if: steps.cache-tsgcn.outputs.cache-hit != 'true'
+#      uses: actions/checkout@v3
+#      with:
+#        repository: ReactionMechanismGenerator/TS-GCN
+#        path: TS-GCN
+#        ref: main
+#        fetch-depth: 1
+
     - name: Cache KinBot
       id: cache-kinbot
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
           path: KinBot
           key: ${{ runner.os }}-kinbot-2.0.6
@@ -122,14 +126,14 @@ jobs:
         fetch-depth: 1
 
     - name: Cache Packages
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       env:
         CACHE_NUMBER: 0
       with:
         path: ~/conda_pkgs_dir
         key:
           ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-${{ hashFiles('environment.yml') }}
-    
+
     - name: Set up miniconda
       uses: conda-incubator/setup-miniconda@v3
       with:
@@ -140,7 +144,7 @@ jobs:
           conda-solver: libmamba
 
     - name: Cache ARC env
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
         path: ${{ env.CONDA }}/envs/arc_env
         key: conda-${{ runner.os }}--${{ runner.arch }}-arcenv-${{ env.CACHE_NUMBER}}
@@ -164,14 +168,14 @@ jobs:
         echo 'export PATH=$PATH:'"$(pwd)" >> ~/.bashrc
         make
         cd ..
-  
+
     - name: Set ARC Path
       run: |
         echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> $GITHUB_ENV
         echo 'export PATH=$PATH:'"$(pwd)" >> $GITHUB_ENV
         echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> ~/.bashrc
         echo 'export PATH=$PATH:'"$(pwd)" >> ~/.bashrc
-        
+
     - name: Install AutoTST
       run: |
         cd AutoTST
@@ -183,18 +187,18 @@ jobs:
         # install pyaml
         conda install -n tst_env -c conda-forge -y pyyaml
         cd ..
-    
-    - name: Install TS-GCN
-      run: |
-        cd TS-GCN
-        echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> $GITHUB_ENV
-        echo 'export PATH=$PATH:'"$(pwd)" >> $GITHUB_ENV
-        echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> ~/.bashrc
-        echo 'export PATH=$PATH:'"$(pwd)" >> ~/.bashrc
-        bash devtools/create_env_cpu.sh
-        conda env update -n ts_gcn -f environment.yml
-        cd ..
-    
+
+#    - name: Install TS-GCN
+#      run: |
+#        cd TS-GCN
+#        echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> $GITHUB_ENV
+#        echo 'export PATH=$PATH:'"$(pwd)" >> $GITHUB_ENV
+#        echo 'export PYTHONPATH=$PYTHONPATH:'"$(pwd)" >> ~/.bashrc
+#        echo 'export PATH=$PATH:'"$(pwd)" >> ~/.bashrc
+#        bash devtools/create_env_cpu.sh
+#        conda env update -n ts_gcn -f environment.yml
+#        cd ..
+
     - name: Install KinBot
       run: |
         cd KinBot
@@ -218,13 +222,13 @@ jobs:
     - name: Install Torch Ani
       run: |
         conda env create -f devtools/tani_environment.yml
-  
+
     - name: Install XTB
       run: |
         conda env create -f devtools/xtb_environment.yml
-    
+
     - name: Test with pytest
-      shell: bash -el {0} 
+      shell: bash -el {0}
       run: |
         source ~/.bashrc
         conda activate arc_env

ARC.py

@@ -21,9 +21,8 @@ def parse_command_line_arguments(command_line_args=None):
         command_line_args: The command line arguments.
 
     Returns:
-        The parsed command-line arguments by key words.
+        The parsed command-line arguments by keywords.
     """
-
     parser = argparse.ArgumentParser(description='Automatic Rate Calculator (ARC)')
     parser.add_argument('file', metavar='FILE', type=str, nargs=1,
                         help='a file describing the job to execute')
@@ -45,7 +44,6 @@ def main():
     """
     The main ARC executable function
     """
-
     args = parse_command_line_arguments()
     input_file = args.file
     project_directory = os.path.abspath(os.path.dirname(args.file))

arc/__init__.py

@@ -12,6 +12,8 @@
 import arc.utils
 
 import arc.job
+import arc.family
+import arc.reaction
 import arc.settings
 import arc.species
 import arc.statmech

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,
@@ -18,11 +17,9 @@
                         read_yaml_file,
                         sum_list_entries,
                         )
+from arc.family.family import get_reaction_family_products
 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.statmech.factory import statmech_factory
 
 if TYPE_CHECKING:
@@ -302,12 +299,10 @@ 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
-    rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=reaction) or list()
+    product_dicts = get_reaction_family_products(rxn=reaction, rmg_family_set=reaction.family if reaction.family else None)
     freqs, normal_modes_disp = parser.parse_normal_mode_displacement(path=job.local_path_to_output_file, raise_error=False)
     if not len(normal_modes_disp):
         return
@@ -332,11 +327,8 @@ def check_normal_mode_displacement(reaction: 'ARCReaction',
                                            tolerance=1.5,
                                            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]
-            if r_label_dict is None:
-                continue
+        for i, product_dict in enumerate(product_dicts):
+            r_label_dict = product_dict['r_label_map']
             expected_breaking_bonds, expected_forming_bonds = reaction.get_expected_changing_bonds(r_label_dict=r_label_dict)
             if expected_breaking_bonds is None or expected_forming_bonds is None:
                 continue
@@ -355,7 +347,7 @@ def check_normal_mode_displacement(reaction: 'ARCReaction',
                                                          'breaking/forming bonds due to a missing RMG template; '
             reaction.ts_species.ts_checks['NMD'] = True
             break
-        if not len(rmg_reactions):
+        if not len(product_dicts):
             # Just check that some bonds break/form, and that this is not a torsional saddle point.
             warning = f'Cannot check normal mode displacement for reaction {reaction} since a corresponding ' \
                       f'RMG template could not be generated'
@@ -535,18 +527,17 @@ def get_rxn_normal_mode_disp_atom_number(rxn_family: Optional[str] = None,
     Returns:
         int: The respective number of atoms.
     """
+    if rxn_family is None and reaction is None:
+        raise ValueError('Either `rxn_family` or `reaction` must be given.')
     default = 3
     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 is not None and 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

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()
@@ -659,7 +646,8 @@ def test_get_rxn_normal_mode_disp_atom_number(self):
             ts.get_rxn_normal_mode_disp_atom_number('family', rms_list=['family'])
         with self.assertRaises(TypeError):
             ts.get_rxn_normal_mode_disp_atom_number('family', rms_list=15.215)
-        self.assertEqual(ts.get_rxn_normal_mode_disp_atom_number(), 3)
+        with self.assertRaises(ValueError):
+            self.assertEqual(ts.get_rxn_normal_mode_disp_atom_number(), 3)
         self.assertEqual(ts.get_rxn_normal_mode_disp_atom_number('default'), 3)
         self.assertEqual(ts.get_rxn_normal_mode_disp_atom_number('intra_H_migration'), 3)
         self.assertEqual(ts.get_rxn_normal_mode_disp_atom_number('intra_H_migration', rms_list=self.rms_list_1), 4)