Merge pull request #1769 from ReactionMechanismGenerator/arkane_libs

Output RMG libraries in Arkane

arkane/input.py

@@ -645,7 +645,7 @@ def load_input_file(path):
             if atom_energies is not None:
                 job.arkane_species.atom_energies = atom_energies
 
-    return job_list, reaction_dict, species_dict, transition_state_dict, network_dict
+    return job_list, reaction_dict, species_dict, transition_state_dict, network_dict, model_chemistry
 
 
 def process_model_chemistry(model_chemistry):

arkane/inputTest.py

@@ -211,7 +211,8 @@ def test_load_input_file(self):
         """Test loading an Arkane input file"""
         path = os.path.join(os.path.dirname(os.path.dirname(rmgpy.__file__)), 'examples', 'arkane', 'networks',
                             'acetyl+O2', 'input.py')
-        job_list, reaction_dict, species_dict, transition_state_dict, network_dict = load_input_file(path)
+        job_list, reaction_dict, species_dict, transition_state_dict, network_dict, model_chemistry \
+            = load_input_file(path)
 
         self.assertEqual(len(job_list), 1)
 
@@ -230,6 +231,8 @@ def test_load_input_file(self):
         self.assertEqual(len(network_dict), 1)
         self.assertTrue('acetyl + O2' in network_dict)
 
+        self.assertEqual(model_chemistry, '')
+
     def test_process_model_chemistry(self):
         """
         Test processing the model chemistry to derive the sp and freq levels

arkane/main.py

@@ -47,6 +47,7 @@
 except ImportError:
     pass
 
+from rmgpy import __version__
 from rmgpy.chemkin import write_elements_section
 from rmgpy.data.thermo import ThermoLibrary
 from rmgpy.data.base import Entry
@@ -57,12 +58,11 @@
 from arkane.explorer import ExplorerJob
 from arkane.input import load_input_file
 from arkane.kinetics import KineticsJob
+from arkane.output import save_thermo_lib, save_kinetics_lib
 from arkane.pdep import PressureDependenceJob
 from arkane.statmech import StatMechJob
 from arkane.thermo import ThermoJob
 
-################################################################################
-
 
 class Arkane(object):
     """
@@ -88,11 +88,12 @@ class Arkane(object):
     :meth:`parse_command_line_arguments()` method before running :meth:`execute()`.
     """
 
-    def __init__(self, input_file=None, output_directory=None, verbose=logging.INFO):
+    def __init__(self, input_file=None, output_directory=None, verbose=logging.INFO, save_rmg_libraries=True):
         self.job_list = []
         self.input_file = input_file
         self.output_directory = output_directory
         self.verbose = verbose
+        self.save_rmg_libraries = save_rmg_libraries
 
     def parse_command_line_arguments(self):
         """
@@ -152,8 +153,8 @@ def load_input_file(self, input_file):
         loaded set of jobs as a list.
         """
         self.input_file = input_file
-        self.job_list, self.reaction_dict, self.species_dict, self.transition_state_dict, self.network_dict = \
-            load_input_file(self.input_file)
+        self.job_list, self.reaction_dict, self.species_dict, self.transition_state_dict, self.network_dict, \
+            self.model_chemistry = load_input_file(self.input_file)
         logging.info('')
         return self.job_list
 
@@ -281,6 +282,25 @@ def execute(self):
         # Print some information to the end of the log
         log_footer()
 
+        if self.save_rmg_libraries:
+            # save RMG thermo and kinetics libraries
+            species, reactions = list(), list()
+            for job in self.job_list:
+                if isinstance(job, ThermoJob) and len(job.species.molecule):
+                    species.append(job.species)
+                elif isinstance(job, KineticsJob) \
+                        and all([len(species.molecule) for species in job.reaction.reactants + job.reaction.products]):
+                    reactions.append(job.reaction)
+                elif isinstance(job, PressureDependenceJob):
+                    for reaction in job.network.path_reactions:
+                        if all([len(species.molecule) for species in reaction.reactants + reaction.products]):
+                            reactions.append(reaction)
+            lib_path = os.path.join(self.output_directory, 'RMG_libraries')
+            model_chemistry = f' at the {self.model_chemistry} level of theory' if self.model_chemistry else ''
+            lib_long_desc = f'Calculated using Arkane v{__version__}{model_chemistry}.'
+            save_thermo_lib(species_list=species, path=lib_path, name='thermo', lib_long_desc=lib_long_desc)
+            save_kinetics_lib(rxn_list=reactions, path=lib_path, name='kinetics', lib_long_desc=lib_long_desc)
+
     def get_libraries(self):
         """Get RMG kinetics and thermo libraries"""
         name = 'kineticsjobs'

arkane/output.py

@@ -32,8 +32,14 @@
 """
 
 import ast
+import logging
+import os
+import shutil
 
-################################################################################
+from rmgpy.data.base import Entry
+from rmgpy.data.kinetics.library import KineticsLibrary
+from rmgpy.data.thermo import ThermoLibrary
+from rmgpy.species import Species
 
 
 class PrettifyVisitor(ast.NodeVisitor):
@@ -185,3 +191,94 @@ def prettify(string, indent=4):
     visitor.visit(node)
     # Return the pretty version of the string
     return visitor.string
+
+
+def get_str_xyz(spc):
+    """
+    Get a string representation of the 3D coordinates from the conformer.
+
+    Args:
+        spc (Species): A Species instance.
+
+    Returns:
+        str: A string representation of the coordinates
+    """
+    if spc.conformer.coordinates is not None:
+        from arkane.common import symbol_by_number
+        xyz_list = list()
+        for number, coord in zip(spc.conformer.number.value_si, spc.conformer.coordinates.value_si):
+            coord_angstroms = coord * 10 ** 10
+            row = f'{symbol_by_number[number]:4}'
+            row += '{0:14.8f}{1:14.8f}{2:14.8f}'.format(*coord_angstroms)
+            xyz_list.append(row)
+        return '\n'.join(xyz_list)
+    else:
+        return None
+
+
+def save_thermo_lib(species_list, path, name, lib_long_desc):
+    """
+    Save an RMG thermo library.
+
+    Args:
+        species_list (list): Entries are Species object instances for which thermo will be saved.
+        path (str): The base folder in which the thermo library will be saved.
+        name (str): The library name.
+        lib_long_desc (str): A multiline string with relevant description.
+    """
+    if species_list:
+        lib_path = os.path.join(path, f'{name}.py')
+        thermo_library = ThermoLibrary(name=name, long_desc=lib_long_desc)
+        for i, spc in enumerate(species_list):
+            if spc.thermo is not None:
+                long_thermo_description = f'\nSpin multiplicity: {spc.conformer.spin_multiplicity}' \
+                                          f'\nExternal symmetry: {spc.molecule[0].symmetry_number}' \
+                                          f'\nOptical isomers: {spc.conformer.optical_isomers}\n'
+                xyz = get_str_xyz(spc)
+                if xyz is not None:
+                    long_thermo_description += f'\nGeometry:\n{xyz}'
+                thermo_library.load_entry(index=i,
+                                          label=spc.label,
+                                          molecule=spc.molecule[0].to_adjacency_list(),
+                                          thermo=spc.thermo,
+                                          shortDesc=spc.thermo.comment,
+                                          longDesc=long_thermo_description)
+            else:
+                logging.warning(f'Species {spc.label} did not contain any thermo data and was omitted from the thermo '
+                                f'library {name}.')
+        thermo_library.save(lib_path)
+
+
+def save_kinetics_lib(rxn_list, path, name, lib_long_desc):
+    """
+    Save an RMG kinetics library.
+
+    Args:
+        rxn_list (list): Entries are Reaction object instances for which kinetics will be saved.
+        path (str): The base folder in which the kinetic library will be saved.
+        name (str): The library name.
+        lib_long_desc (str): A multiline string with relevant description.
+    """
+    entries = dict()
+    if rxn_list:
+        for i, rxn in enumerate(rxn_list):
+            if rxn.kinetics is not None:
+                entry = Entry(
+                    index=i,
+                    item=rxn,
+                    data=rxn.kinetics,
+                    label=rxn.label)
+                entries[i+1] = entry
+            else:
+                logging.warning(f'Reaction {rxn.label} did not contain any kinetic data and was omitted from the '
+                                f'kinetics library.')
+        kinetics_library = KineticsLibrary(name=name, long_desc=lib_long_desc, auto_generated=True)
+        kinetics_library.entries = entries
+        if os.path.exists(path):
+            shutil.rmtree(path)
+        try:
+            os.makedirs(path)
+        except OSError:
+            pass
+        kinetics_library.save(os.path.join(path, 'reactions.py'))
+        kinetics_library.save_dictionary(os.path.join(path, 'dictionary.txt'))

arkane/outputTest.py

@@ -39,8 +39,11 @@
 from nose.plugins.attrib import attr
 
 import rmgpy
+
+from arkane.gaussian import GaussianLog
 from arkane.main import Arkane
-from arkane.output import prettify
+from arkane.output import prettify, get_str_xyz
+from rmgpy.species import Species
 
 
 @attr('functional')
@@ -129,6 +132,23 @@ def test_prettify(self):
 )"""
         self.assertEqual(prettify(input_str), expected_output)
 
+    def test_get_str_xyz(self):
+        """Test generating an xyz string from the species.conformer object"""
+        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'ethylene_G3.log'))
+        conformer = log.load_conformer()[0]
+        coords, number, mass = log.load_geometry()
+        conformer.coordinates, conformer.number, conformer.mass = (coords, "angstroms"), number, (mass, "amu")
+        spc1 = Species(smiles='C=C')
+        spc1.conformer = conformer
+        xyz_str = get_str_xyz(spc1)
+        expected_xyz_str = """C       0.00545100    0.00000000    0.00339700
+H       0.00118700    0.00000000    1.08823200
+H       0.97742900    0.00000000   -0.47841600
+C      -1.12745800    0.00000000   -0.70256500
+H      -1.12319800    0.00000000   -1.78740100
+H      -2.09943900    0.00000000   -0.22075700"""
+        self.assertEqual(xyz_str, expected_xyz_str)
+
 
 if __name__ == '__main__':
     unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))