Liquid reactor concentration accepts Quantity init

It's very inconvenient to make a round trip read/write of RMG input
files when the units of the concentration aren't saved in the reactor
system. So this updates the liquid reactor (and input tests) to accept
Quantity types. It still accepts floats, but RMG's input module now
passes it in as a quantity.

This also fixes an error where someone forgot to use deepcopy and so the
Temperature was getting added to the dictionary of initial conditions
instead of just to sensitive conditions

Author: sevyharris

rmgpy/rmg/input.py

@@ -984,7 +984,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,

rmgpy/solver/liquid.pyx

@@ -73,7 +73,10 @@ cdef class LiquidReactor(ReactionSystem):
             self.Trange = [Quantity(t) for t in T]
 
         self.P = Quantity(100000., 'kPa')  # Arbitrary high pressure (1000 Bar) to get reactions in the high-pressure limit!
-        self.initial_concentrations = initial_concentrations  # should be passed in SI
+        for spec, conc in initial_concentrations.items():  # convert to Quantity if given as float (SI mol/m^3)
+            if isinstance(initial_concentrations[spec], float):
+                initial_concentrations[spec] = ScalarQuantity(initial_concentrations[spec], 'mol/m^3')
+        self.initial_concentrations = initial_concentrations  # should be passed in as Quantity to preserve units
         self.V = 0  # will be set from initial_concentrations in initialize_model
         self.constant_volume = True
         self.viscosity = 0  # in Pa*s
@@ -91,8 +94,6 @@ cdef class LiquidReactor(ReactionSystem):
         """
         initial_concentrations = {}
         for label, moleFrac in self.initial_concentrations.items():
-            if label == 'T':
-                continue
             initial_concentrations[species_dict[label]] = moleFrac
         self.initial_concentrations = initial_concentrations
 
@@ -202,7 +203,7 @@ cdef class LiquidReactor(ReactionSystem):
 
         for spec, conc in self.initial_concentrations.items():
             i = self.get_species_index(spec)
-            self.core_species_concentrations[i] = conc
+            self.core_species_concentrations[i] = conc.value_si
 
         V = 1.0 / np.sum(self.core_species_concentrations)
         self.V = V

test/rmgpy/rmg/inputTest.py

@@ -280,10 +280,18 @@ def test_liquid_reactor_concentrations(self):
         global rmg
         reactor = rmg.reaction_systems[0]
 
-        # Values get converted to default SI units, mol/m^3
-        assert reactor.initial_concentrations["A"] == 300
-        assert reactor.initial_concentrations["B"] == 200
-        assert reactor.initial_concentrations["C"] == 100
+        # Values are saved as Quantity to preserve units, but get changed to SI units, mol/m^3 in the set_initial_conditions function
+        assert reactor.initial_concentrations["A"].value_si == 300
+        assert reactor.initial_concentrations["B"].value_si == 200
+        assert reactor.initial_concentrations["C"].value_si == 100
+
+        rmg.reaction_systems[0].initialize_model(list(reactor.initial_concentrations.keys()), [], [], [])
+        rmg.reaction_systems[0].set_initial_conditions()
+        assert rmg.reaction_systems[0].core_species_concentrations[0] == 300
+        assert rmg.reaction_systems[0].core_species_concentrations[1] == 200
+        assert rmg.reaction_systems[0].core_species_concentrations[2] == 100
+        # Test that the units are saved internally in the reactor as mol/m^3
+
 
     def test_surface_reactor_mole_fractions(self):
         """Test that SurfaceReactor mole fractions are set properly"""