Merge pull request #463 from ReactionMechanismGenerator/remove_1st_atom

Added a function to remove the 1st atom in a Z matrix

arc/common.py

@@ -196,7 +196,7 @@ def initialize_log(log_file: str,
         project_directory (str, optional): The path to the project directory.
         verbose (int, optional): Specify the amount of log text seen.
     """
-    # backup and delete an existing log file if needed
+    # Backup and delete an existing log file if needed.
     if project_directory is not None and os.path.isfile(log_file):
         if not os.path.isdir(os.path.join(project_directory, 'log_and_restart_archive')):
             os.mkdir(os.path.join(project_directory, 'log_and_restart_archive'))
@@ -208,35 +208,35 @@ def initialize_log(log_file: str,
     logger.setLevel(verbose)
     logger.propagate = False
 
-    # Use custom level names for cleaner log output
+    # Use custom level names for cleaner log output.
     logging.addLevelName(logging.CRITICAL, 'Critical: ')
     logging.addLevelName(logging.ERROR, 'Error: ')
     logging.addLevelName(logging.WARNING, 'Warning: ')
     logging.addLevelName(logging.INFO, '')
     logging.addLevelName(logging.DEBUG, '')
     logging.addLevelName(0, '')
 
-    # Create formatter and add to handlers
+    # Create formatter and add to handlers.
     formatter = logging.Formatter('%(levelname)s%(message)s')
 
-    # Remove old handlers before adding ours
+    # Remove old handlers before adding ours.
     while logger.handlers:
         logger.removeHandler(logger.handlers[0])
 
-    # Create console handler; send everything to stdout rather than stderr
+    # Create console handler; send everything to stdout rather than stderr.
     ch = logging.StreamHandler(sys.stdout)
     ch.setLevel(verbose)
     ch.setFormatter(formatter)
     logger.addHandler(ch)
 
-    # Create file handler
+    # Create file handler.
     fh = logging.FileHandler(filename=log_file)
     fh.setLevel(verbose)
     fh.setFormatter(formatter)
     logger.addHandler(fh)
     log_header(project=project)
 
-    # ignore Paramiko, cclib, and matplotlib warnings:
+    # Ignore Paramiko, cclib, and matplotlib warnings:
     warnings.filterwarnings(action='ignore', module='.*paramiko.*')
     warnings.filterwarnings(action='ignore', module='.*cclib.*')
     warnings.filterwarnings(action='ignore', module='.*matplotlib.*')
@@ -258,7 +258,7 @@ def log_header(project: str,
 
     Args:
         project (str): The ARC project name to be logged in the header.
-        level: The desired logging level.
+        level (int, optional): The desired logging level.
     """
     logger.log(level, f'ARC execution initiated on {time.asctime()}')
     logger.log(level, '')
@@ -274,7 +274,7 @@ def log_header(project: str,
 
     paths_dict = {'ARC': ARC_PATH, 'RMG-Py': RMG_PATH, 'RMG-database': RMG_DATABASE_PATH}
     for repo, path in paths_dict.items():
-        # Extract HEAD git commit
+        # Extract HEAD git commit.
         head, date = get_git_commit(path)
         branch_name = get_git_branch(path)
         if head != '' and date != '':
@@ -296,7 +296,7 @@ def log_footer(execution_time: str,
 
     Args:
         execution_time (str): The overall execution time for ARC.
-        level: The desired logging level.
+        level (int, optional): The desired logging level.
     """
     logger.log(level, '')
     logger.log(level, f'Total execution time: {execution_time}')

arc/rmgdb.py

@@ -83,8 +83,8 @@ def load_rmg_database(rmgdb: RMGDatabase,
         include_nist (bool, optional): Whether to include the NIST kinetics libraries,
                                       ``True`` to include, default is ``False``
     """
-    thermo_libraries = thermo_libraries if thermo_libraries is not None else list()
-    reaction_libraries = reaction_libraries if reaction_libraries is not None else list()
+    thermo_libraries = thermo_libraries or list()
+    reaction_libraries = reaction_libraries or list()
     if isinstance(thermo_libraries, str):
         thermo_libraries.replace(' ', '')
         thermo_libraries = [lib for lib in thermo_libraries.split(',')]
@@ -117,22 +117,16 @@ def load_rmg_database(rmgdb: RMGDatabase,
         thermo_libraries = thermo_priority + thermo_libraries
 
     if not reaction_libraries:
-        kinetics_path = os.path.join(db_path, 'kinetics', 'libraries')
-        # Avoid reading .DS_store files for compatibility with Mac OS
-        reaction_libraries = [library for library in os.listdir(kinetics_path) if not library.startswith('.')]
-        indices_to_pop = list()
-        second_level_libraries = list()
-        for i, library in enumerate(reaction_libraries):
-            if not os.path.isfile(os.path.join(kinetics_path, library, 'reactions.py')) and library != 'Surface':
-                indices_to_pop.append(i)
-                # Avoid reading .DS_store files for compatibility with Mac OS (`if not second_level.startswith('.')`)
-                second_level_libraries.extend([library + '/' + second_level
-                                               for second_level in os.listdir(os.path.join(kinetics_path, library))
-                                               if not second_level.startswith('.')])
-        for i in reversed(range(len(reaction_libraries))):  # pop starting from the end, so other indices won't change
-            if i in indices_to_pop or reaction_libraries[i] == 'Surface':
-                reaction_libraries.pop(i)
-        reaction_libraries.extend(second_level_libraries)
+        libraries_path = os.path.join(db_path, 'kinetics', 'libraries')
+        for path, folders, _ in os.walk(libraries_path):
+            if 'Surface' in path:
+                continue
+            for folder in folders:
+                if os.path.isfile(os.path.join(path, folder, 'reactions.py')):
+                    prefix = path.split(libraries_path)[1]
+                    library = os.path.join(prefix, folder)
+                    library = library[1:] if library[0] == '/' else library
+                    reaction_libraries.append(library)
 
     if not load_kinetic_libs:
         reaction_libraries = list()
@@ -312,7 +306,7 @@ def determine_rmg_kinetics(rmgdb: RMGDatabase,
     Returns: List[Reaction]
         All matching RMG reactions (both libraries and families) with a populated ``.kinetics`` attribute.
     """
-    reaction = reaction.copy()  # use a copy to avoid changing atom order in the molecules by RMG
+    reaction = reaction.copy()
     rmg_reactions = list()
     # Libraries:
     for library in rmgdb.kinetics.libraries.values():

arc/species/conformersTest.py

@@ -138,11 +138,11 @@ def test_generate_conformers_with_openbabel(self):
         O          1.39175        1.59510        0.11494""")
         spc = ARCSpecies(label='CCO[O]', smiles='CCO[O]', xyz=xyz)
         lowest_confs = conformers.generate_conformers(mol_list=spc.mol_list, label=spc.label,
-                                                        charge=spc.charge, multiplicity=spc.multiplicity,
-                                                        force_field='MMFF94s', print_logs=False, diastereomers=None,
-                                                        n_confs=1, return_all_conformers=False)
+                                                      charge=spc.charge, multiplicity=spc.multiplicity,
+                                                      force_field='MMFF94s', print_logs=False, diastereomers=None,
+                                                      n_confs=1, return_all_conformers=False)
         self.assertEqual(len(lowest_confs), 1)
-        self.assertAlmostEqual(lowest_confs[0]['FF energy'], -0.7418115656748858)
+        self.assertAlmostEqual(lowest_confs[0]['FF energy'], -0.74181, 2)
         expected_xyz = {'symbols': ('C', 'C', 'H', 'H', 'H', 'O', 'H', 'H', 'O'),
                         'isotopes': (12, 12, 1, 1, 1, 16, 1, 1, 16),
                         'coords': ((-1.06401, 0.15134, -0.02907),
@@ -576,8 +576,8 @@ def test_openbabel_force_field_on_rdkit_conformers(self):
                                     (0.59336, 1.16044, 0.53412),
                                     (2.51642, -0.38757, 0.25036))}]
         self.assertEqual(len(energies), 2)
-        self.assertAlmostEqual(energies[0], 3.4958104235174376)
-        self.assertAlmostEqual(energies[1], -0.7028477774918569)
+        self.assertAlmostEqual(energies[0], 3.49581, 5)
+        self.assertAlmostEqual(energies[1], -0.702848, 5)
         # Only symbols instead of the coordinate values are compared.
         # This is due to the unknown behavior of OpenBabel optimization function.
         # With the same iteration number and same initial xyz, the optimized xyzs can

arc/species/speciesTest.py

@@ -556,7 +556,8 @@ def test_as_dict(self):
                          'is_ts': False,
                          'bond_corrections': {'C-H': 3, 'C-N': 1, 'H-N': 2}}
         self.assertEqual(spc_dict, expected_dict)
-        self.assertTrue(spc_dict['mol']['atoms'][0]['id'] < -50)
+        self.assertEqual(len(set([spc_dict['mol']['atoms'][i]['id'] for i in range(len(spc_dict['mol']['atoms']))])),
+                         len(spc_dict['mol']['atoms']))
         self.assertEqual(spc_dict['mol']['atoms'][0]['radical_electrons'], 0)
         self.assertEqual(spc_dict['mol']['atoms'][0]['charge'], 0)
         self.assertEqual(spc_dict['mol']['atoms'][0]['lone_pairs'], 0)