implement writing of input files with liquidSurfaceReactor

Author: sevyharris

rmgpy/rmg/input.py

@@ -58,6 +58,7 @@
 from rmgpy.solver.mbSampled import MBSampledReactor
 from rmgpy.solver.simple import SimpleReactor
 from rmgpy.solver.surface import SurfaceReactor
+from rmgpy.solver.base import ReactionSystem
 from rmgpy.solver.termination import (
     TerminationConversion,
     TerminationRateRatio,
@@ -656,7 +657,10 @@ def liquid_cat_reactor(temperature,
                    constantSpecies=[]):
     for spec, conc in initialConcentrations.items():
         if not isinstance(conc, list):
-            initialConcentrations[spec] = Quantity(conc)
+            concentration = Quantity(conc)
+            # check the dimensions are ok
+            # convert to mol/m^3 (or something numerically nice? or must it be SI)
+            initialConcentrations[spec] = concentration.value_si
         else:
             if len(conc) != 2:
                 raise InputError("Concentration values must either be in the form of (number,units) or a list with 2 "
@@ -712,14 +716,14 @@ def liquid_cat_reactor(temperature,
 
     initialCondLiq = dict()
     V = 1.0
-    A = V * Quantity(surfaceVolumeRatio).value_si
-    for key, conc in initialConcentrations.items():
-        initialCondLiq[key] = conc.value_si * V
+    A = V*Quantity(surfaceVolumeRatio).value_si
+    for key,item in initialConcentrations.items():
+        initialCondLiq[key] = item*V
     initialCondLiq["T"] = T
     initialCondLiq["V"] = V
     initialCondSurf = dict()
-    for key, surf_cov in initialSurfaceCoverages.items():
-        initialCondSurf[key] = surf_cov * rmg.surface_site_density.value_si * A
+    for key,item in initialSurfaceCoverages.items():
+        initialCondSurf[key] = item*rmg.surface_site_density.value_si*A
     initialCondSurf["T"] = T
     initialCondSurf["A"] = A
     initialCondSurf["d"] = 0.0
@@ -768,7 +772,8 @@ def constant_T_V_liquid_reactor(temperature,
 
     for spec, conc in initialConcentrations.items():
         if not isinstance(conc, list):
-            initialConcentrations[spec] = Quantity(conc)
+            concentration = Quantity(conc)
+            initialConcentrations[spec] = concentration.value_si
         else:
             raise InputError("Condition ranges not supported for this reaction type")
             if len(conc) != 2:
@@ -876,16 +881,16 @@ def constant_T_V_liquid_reactor(temperature,
     ############################################### process inputs ##############################################
 
     initial_conditions = dict()
-    for key, conc in initialConcentrations.items():
-        initial_conditions[key] = conc.value_si * V
+    for key, item in initialConcentrations.items():
+        initial_conditions[key] = item*V
     initial_conditions["T"] = T
     initial_conditions["V"] = V
 
     inlet_conditions = dict()
     if inletConcentrations:
         total_molar_flow_rate = 0
-        for key, inlet_conc in inletConcentrations.items():
-            inlet_conditions[key] = inlet_conc.value_si * inlet_volumetric_flow_rate
+        for key, item in inletConcentrations.items():
+            inlet_conditions[key] = item*inlet_volumetric_flow_rate
             total_molar_flow_rate += inlet_conditions[key]
         for key, item in inlet_conditions.items():
             inlet_conditions[key] = item/total_molar_flow_rate #molar fraction for each species
@@ -942,7 +947,10 @@ def liquid_reactor(temperature,
 
     for spec, conc in initialConcentrations.items():
         if not isinstance(conc, list):
-            initialConcentrations[spec] = Quantity(conc)
+            concentration = Quantity(conc)
+            # check the dimensions are ok
+            # convert to mol/m^3 (or something numerically nice? or must it be SI)
+            initialConcentrations[spec] = concentration.value_si
         else:
             if len(conc) != 2:
                 raise InputError("Concentration values must either be in the form of (number,units) or a list with 2 "
@@ -983,7 +991,7 @@ def liquid_reactor(temperature,
     if sensitivityConcentrations is None or sensitivityTemperature is None:
         sens_conditions = None
     else:
-        sens_conditions = sensitivityConcentrations
+        sens_conditions = deepcopy(sensitivityConcentrations)
         sens_conditions['T'] = Quantity(sensitivityTemperature).value_si
 
     system = LiquidReactor(T, initialConcentrations, nSims, termination, sensitive_species, sensitivityThreshold,
@@ -1749,7 +1757,29 @@ def format_initial_mole_fractions(system):
 
     # Reaction systems
     for system in rmg.reaction_systems:
-        if rmg.solvent:
+        if isinstance(system, ConstantTLiquidSurfaceReactor):
+            f.write('liquidSurfaceReactor(\n')
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    initialConcentrations={\n')
+            for spcs, conc in system.initial_conditions['liquid'].items():
+                if spcs in ['T', 'V']:
+                    continue
+                f.write('        "{0!s}": ({1:g},"{2!s}"),\n'.format(spcs, conc, 'mol/m^3'))
+            f.write('    initialSurfaceCoverages={\n')
+            for spcs, conc_mols in system.initial_conditions['surface'].items():
+                if spcs in ['T', 'A', 'd']:
+                    continue
+                # surf conc here is in mols, need to convert back into unitless coverage fraction
+                coverage = conc_mols / (rmg.surface_site_density.value_si * system.initial_conditions['surface']['A'])
+                f.write('        "{0!s}": ({1:g}),\n'.format(spcs, coverage))
+            f.write('    },\n')
+            
+            # write the list of constant species
+            f.write(f'    constantSpecies = {system.const_spc_names},\n')
+
+            # write the surface Volume ratio, where ratio = A/V and A was originally constructed by assuming V=1 m^3
+            f.write('    surfaceVolumeRatio = ({0:g}, "{1!s}"),\n'.format(system.initial_conditions['surface']['A'], 'm^-1'))
+        elif isinstance(system, LiquidReactor):
             f.write('liquidReactor(\n')
             f.write('    temperature = ' + format_temperature(system) + '\n')
             f.write('    initialConcentrations={\n')
@@ -1778,8 +1808,25 @@ def format_initial_mole_fractions(system):
         f.write('    },\n')
 
         # Termination criteria
+        if isinstance(system, ReactionSystem):
+            terminations = system.termination
+        elif isinstance(system, Reactor):  # RMS reactor terminations need to be converted back
+            terminations = []
+            for term in system.terminations:
+                if hasattr(term, 'time'):
+                    terminations.append(TerminationTime(time=(term.time, 's')))
+                elif hasattr(term, 'ratio'):
+                    terminations.append(TerminationRateRatio(ratio=term.ratio))
+                elif isinstance(term, tuple):
+                    species, conversion = term
+                    terminations.append(TerminationConversion(spec=species, conv=conversion))
+                else:
+                    raise NotImplementedError('Termination criterion of type {0} is not currently supported for RMS reactors. Please convert this criterion to a time-based criterion or remove it from the input file.'.format(type(term)))
+        else:
+            raise NotImplementedError('Termination criteria for reaction system of type {0} not supported'.format(type(system)))
+
         conversions = ''
-        for term in system.termination:
+        for term in terminations:
             if isinstance(term, TerminationTime):
                 f.write('    terminationTime = ({0:g},"{1!s}"),\n'.format(term.time.value, term.time.units))
             elif isinstance(term, TerminationRateRatio):