Merge pull request #2532 to fix rmgpy.rmg.input.save_input_file

Fix writing of RMG input files by rmgpy.rmg.input.save_input_file
A few fixes, so that scripts can use RMG to (load, modify, and) save input files.
(At least for many more situations than they formerly could, even if we're not at 100% yet.)

Author: rwest

rmgpy/rmg/input.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+
 
 ###############################################################################
 #                                                                             #
@@ -44,7 +44,7 @@
 from rmgpy.molecule import Molecule
 from rmgpy.molecule.fragment import Fragment
 from rmgpy.molecule.group import Group
-from rmgpy.quantity import Energy, Quantity, RateCoefficient, SurfaceConcentration
+from rmgpy.quantity import Energy, Quantity, RateCoefficient, SurfaceConcentration, Concentration
 from rmgpy.rmg.model import CoreEdgeReactionModel
 from rmgpy.rmg.reactionmechanismsimulator_reactors import (
     ConstantTLiquidSurfaceReactor,
@@ -1643,6 +1643,8 @@ def read_thermo_input_file(path, rmg0):
     rmg.reaction_model = CoreEdgeReactionModel()
     rmg.initial_species = []
     rmg.reaction_systems = []
+    if rmg.output_directory is None:
+        rmg.output_directory = os.path.dirname(full_path)
     species_dict = {}
 
     global_context = {'__builtins__': None}
@@ -1699,12 +1701,16 @@ def save_input_file(path, rmg):
     f.write('    kineticsEstimator = {0!r},\n'.format(rmg.kinetics_estimator))
     f.write(')\n\n')
 
-    if rmg.surfaceSiteDenisty or rmg.binding_energies:
+    if rmg.surface_site_density or rmg.binding_energies:
         f.write('catalystProperties(\n')
-        if rmg.surfaceSiteDenisty:
-            f.write('    surface_site_density = {0!r},'.format(rmg.surface_site_density))
+        if rmg.surface_site_density:
+            f.write('    surfaceSiteDensity = ({0:g}, "{1!s}"),'.format(rmg.surface_site_density.value, rmg.surface_site_density.units) + '\n')
         if rmg.binding_energies:
-            f.write('    binding_energies = {0!r},'.format(rmg.binding_energies))
+            f.write('    bindingEnergies = {\n')
+            for spc, be in rmg.binding_energies.items():
+                f.write('        "{0!s}": ({1:g}, "{2!s}"),\n'.format(spc, be.value, be.units))
+            f.write('    },\n')
+
         f.write(')\n\n')
 
     # Species
@@ -1718,33 +1724,72 @@ def save_input_file(path, rmg):
         f.write('"""),\n')
         f.write(')\n\n')
 
+    def format_temperature(system):
+        """Get temperature string format for reaction system, whether single value or range"""
+        if system.T is not None:
+            return '({0:g},"{1!s}")'.format(system.T.value, system.T.units)
+        
+        return f'[({system.Trange[0].value:g}, "{system.Trange[0].units}"), ({system.Trange[1].value:g}, "{system.Trange[1].units}")],'
+
+    def format_pressure(system):
+        """Get pressure string format for reaction system, whether single value or range"""
+        if system.P is not None:
+            return '({0:g},"{1!s}")'.format(system.P.value, system.P.units)
+        
+        return f'[({system.Prange[0].value:g}, "{system.Prange[0].units}"), ({system.Prange[1].value:g}, "{system.Prange[1].units}")],'
+
+    def format_initial_mole_fractions(system):
+        """Get initial mole fractions string format for reaction system"""
+        mole_fractions = ''
+        for spcs, molfrac in system.initial_mole_fractions.items():
+            if isinstance(molfrac, list):
+                mole_fractions += '        "{0!s}": [{1:g}, {2:g}],\n'.format(spcs.label, molfrac[0], molfrac[1])
+            else:
+                mole_fractions += '        "{0!s}": {1:g},\n'.format(spcs.label, molfrac)
+        return mole_fractions
+
+
     # Reaction systems
     for system in rmg.reaction_systems:
         if rmg.solvent:
             f.write('liquidReactor(\n')
-            f.write('    temperature = ({0:g},"{1!s}"),\n'.format(system.T.value, system.T.units))
+            f.write('    temperature = ' + format_temperature(system) + '\n')
             f.write('    initialConcentrations={\n')
             for spcs, conc in system.initial_concentrations.items():
+                # conc may have been converted to SI, so we need to convert back
+                if type(conc) == float:
+                    conc = Quantity(conc, Concentration.units)
                 f.write('        "{0!s}": ({1:g},"{2!s}"),\n'.format(spcs.label, conc.value, conc.units))
+        elif isinstance(system, SurfaceReactor):
+            f.write('surfaceReactor(\n')
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    initialPressure = ({0:g},"{1!s}"),\n'.format(system.P_initial.value, system.P_initial.units))
+            f.write('    initialGasMoleFractions={\n')
+            for spcs, molfrac in system.initial_gas_mole_fractions.items():
+                f.write('        "{0!s}": {1:g},\n'.format(spcs.label, molfrac))
+            f.write('    },\n')
+            f.write('    initialSurfaceCoverages={\n')
+            for spcs, cov in system.initial_surface_coverages.items():
+                f.write('        "{0!s}": {1:g},\n'.format(spcs.label, cov))
         else:
             f.write('simpleReactor(\n')
-            f.write('    temperature = ({0:g},"{1!s}"),\n'.format(system.T.value, system.T.units))
-            # Convert the pressure from SI pascal units to bar here
-            # Do something more fancy later for converting to user's desired units for both T and P..
-            f.write('    pressure = ({0:g},"{1!s}"),\n'.format(system.P.value, system.P.units))
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    pressure = ' + format_pressure(system) + '\n')
             f.write('    initialMoleFractions={\n')
-            for spcs, molfrac in system.initial_mole_fractions.items():
-                f.write('        "{0!s}": {1:g},\n'.format(spcs.label, molfrac))
+            f.write(format_initial_mole_fractions(system))
         f.write('    },\n')
 
         # Termination criteria
         conversions = ''
         for term in system.termination:
             if isinstance(term, TerminationTime):
                 f.write('    terminationTime = ({0:g},"{1!s}"),\n'.format(term.time.value, term.time.units))
-
-            else:
+            elif isinstance(term, TerminationRateRatio):
+                f.write('    terminationRateRatio = {0:g},\n'.format(term.ratio))
+            elif isinstance(term, TerminationConversion):
                 conversions += '        "{0:s}": {1:g},\n'.format(term.species.label, term.conversion)
+            else:
+                raise NotImplementedError('Termination criteria of type {0} not supported'.format(type(term)))
         if conversions:
             f.write('    terminationConversion = {\n')
             f.write(conversions)
@@ -1787,9 +1832,9 @@ def save_input_file(path, rmg):
     if rmg.pressure_dependence:
         f.write('pressureDependence(\n')
         f.write('    method = {0!r},\n'.format(rmg.pressure_dependence.method))
-        f.write('    maximumGrainSize = ({0:g},"{1!s}"),\n'.format(rmg.pressure_dependence.grain_size.value,
-                                                                   rmg.pressure_dependence.grain_size.units))
-        f.write('    minimumNumberOfGrains = {0},\n'.format(rmg.pressure_dependence.grain_count))
+        f.write('    maximumGrainSize = ({0:g},"{1!s}"),\n'.format(rmg.pressure_dependence.maximum_grain_size.value,
+                                                                   rmg.pressure_dependence.maximum_grain_size.units))
+        f.write('    minimumNumberOfGrains = {0},\n'.format(rmg.pressure_dependence.minimum_grain_count))
         f.write('    temperatures = ({0:g},{1:g},"{2!s}",{3:d}),\n'.format(
             rmg.pressure_dependence.Tmin.value,
             rmg.pressure_dependence.Tmax.value,