diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index f4780e2843..d8dd5807bc 100644
--- a/rmgpy/rmg/input.py
+++ b/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,
@@ -466,7 +467,7 @@ def constant_V_ideal_gas_reactor(temperature,
                 raise InputError('Initial mole fraction range out of order: {0}'.format(key))
 
     if not isinstance(temperature, list):
-        T = Quantity(temperature).value_si
+        T = Quantity(temperature)
     else:
         raise InputError("Condition ranges not supported for this reaction type")
         if len(temperature) != 2:
@@ -514,11 +515,11 @@ def constant_V_ideal_gas_reactor(temperature,
         termination.append(TerminationRateRatio(terminationRateRatio))
     if len(termination) == 0:
         raise InputError('No termination conditions specified for reaction system #{0}.'.format(len(rmg.reaction_systems) + 2))
-    
+
     initial_cond = initialMoleFractions
-    initial_cond["T"] = T 
+    initial_cond["T"] = T.value_si
     initial_cond["P"] = P
-    system = ConstantVIdealGasReactor(rmg.reaction_model.core.phase_system,rmg.reaction_model.edge.phase_system,initial_cond,termination)
+    system = ConstantVIdealGasReactor(rmg.reaction_model.core.phase_system, rmg.reaction_model.edge.phase_system, initial_cond,termination)
     system.T = Quantity(T)
     system.P = Quantity(P)
     system.Trange = None
@@ -690,7 +691,7 @@ def liquid_cat_reactor(temperature,
                 raise InputError('Species {0} not found in the input file'.format(const_spc))
 
     if not isinstance(temperature, list):
-        T = Quantity(temperature).value_si
+        T = Quantity(temperature)
     else:
         raise InputError("Condition ranges not supported for this reaction type")
         if len(temperature) != 2:
@@ -720,12 +721,12 @@ def liquid_cat_reactor(temperature,
     A = V*Quantity(surfaceVolumeRatio).value_si
     for key,item in initialConcentrations.items():
         initialCondLiq[key] = item*V
-    initialCondLiq["T"] = T
+    initialCondLiq["T"] = T.value_si
     initialCondLiq["V"] = V
     initialCondSurf = dict()
     for key,item in initialSurfaceCoverages.items():
         initialCondSurf[key] = item*rmg.surface_site_density.value_si*A
-    initialCondSurf["T"] = T
+    initialCondSurf["T"] = T.value_si
     initialCondSurf["A"] = A
     initialCondSurf["d"] = 0.0
     if surfPotential:
@@ -735,12 +736,12 @@ def liquid_cat_reactor(temperature,
     if distance:
         initialCondLiq["d"] = Quantity(distance).value_si
     if viscosity:
-        initialCondLiq["mu"] = Quantity(distance).value_si
+        initialCondLiq["mu"] = Quantity(viscosity).value_si
     system = ConstantTLiquidSurfaceReactor(rmg.reaction_model.core.phase_system,
                                            rmg.reaction_model.edge.phase_system,
                                            {"liquid":initialCondLiq,"surface":initialCondSurf},
                                            termination,constantSpecies)
-    system.T = Quantity(T)
+    system.T = T
     system.Trange = None
     system.sensitive_species = []
     rmg.reaction_systems.append(system)
@@ -763,7 +764,7 @@ def constant_T_V_liquid_reactor(temperature,
     ################################################# check input ########################################################
 
     if not isinstance(temperature, list):
-        T = Quantity(temperature).value_si
+        T = Quantity(temperature)
     else:
         raise InputError("Condition ranges not supported for this reaction type")
         if len(temperature) != 2:
@@ -884,7 +885,7 @@ def constant_T_V_liquid_reactor(temperature,
     initial_conditions = dict()
     for key, item in initialConcentrations.items():
         initial_conditions[key] = item*V
-    initial_conditions["T"] = T
+    initial_conditions["T"] = T.value_si
     initial_conditions["V"] = V
 
     inlet_conditions = dict()
@@ -919,7 +920,7 @@ def constant_T_V_liquid_reactor(temperature,
                                            outlet_conditions,
                                            evap_cond_conditions,
                                            )
-    system.T = Quantity(T)
+    system.T = T
     system.Trange = None
     system.sensitive_species = []
     rmg.reaction_systems.append(system)
@@ -992,7 +993,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,
@@ -1412,6 +1413,7 @@ def options(name='Seed', generateSeedEachIteration=True, saveSeedToDatabase=Fals
 def generated_species_constraints(**kwargs):
     valid_constraints = [
         'allowed',
+        'explicitlyAllowedMolecules',
         'maximumCarbonAtoms',
         'maximumOxygenAtoms',
         'maximumNitrogenAtoms',
@@ -1575,9 +1577,9 @@ def read_input_file(path, rmg0):
         'adjacencyListGroup': adjacency_list_group,
         'react': react,
         'simpleReactor': simple_reactor,
-        'constantVIdealGasReactor' : constant_V_ideal_gas_reactor,
-        'constantTPIdealGasReactor' : constant_TP_ideal_gas_reactor,
-        'liquidSurfaceReactor' : liquid_cat_reactor,
+        'constantVIdealGasReactor': constant_V_ideal_gas_reactor,
+        'constantTPIdealGasReactor': constant_TP_ideal_gas_reactor,
+        'liquidSurfaceReactor': liquid_cat_reactor,
         'constantTVLiquidReactor': constant_T_V_liquid_reactor,
         'liquidReactor': liquid_reactor,
         'surfaceReactor': surface_reactor,
@@ -1621,6 +1623,8 @@ def read_input_file(path, rmg0):
         if not isinstance(reaction_system, Reactor):
             reaction_system.convert_initial_keys_to_species_objects(species_dict)
 
+    if rmg.output_directory is None:
+        rmg.output_directory = os.path.dirname(full_path)
     if rmg.quantum_mechanics:
         rmg.quantum_mechanics.set_default_output_directory(rmg.output_directory)
         rmg.quantum_mechanics.initialize()
@@ -1705,7 +1709,7 @@ def save_input_file(path, rmg):
     f.write('    thermoLibraries = {0!r},\n'.format(rmg.thermo_libraries))
     f.write('    reactionLibraries = {0!r},\n'.format(rmg.reaction_libraries))
     f.write('    seedMechanisms = {0!r},\n'.format(rmg.seed_mechanisms))
-    f.write('    kinetics_depositories = {0!r},\n'.format(rmg.kinetics_depositories))
+    f.write('    kineticsDepositories = {0!r},\n'.format(rmg.kinetics_depositories))
     f.write('    kineticsFamilies = {0!r},\n'.format(rmg.kinetics_families))
     f.write('    kineticsEstimator = {0!r},\n'.format(rmg.kinetics_estimator))
     f.write(')\n\n')
@@ -1736,15 +1740,15 @@ def save_input_file(path, rmg):
     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 '({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 '({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):
@@ -1757,10 +1761,32 @@ def format_initial_mole_fractions(system):
                 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:
+        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('    },\n')
+            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')
@@ -1769,6 +1795,7 @@ def format_initial_mole_fractions(system):
                 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))
+            f.write('    },\n')
         elif isinstance(system, SurfaceReactor):
             f.write('surfaceReactor(\n')
             f.write('    temperature = ' + format_temperature(system) + '\n')
@@ -1780,17 +1807,64 @@ def format_initial_mole_fractions(system):
             f.write('    initialSurfaceCoverages={\n')
             for spcs, cov in system.initial_surface_coverages.items():
                 f.write('        "{0!s}": {1:g},\n'.format(spcs.label, cov))
+            f.write('    },\n')
+            f.write('    surfaceVolumeRatio = ({0:g}, "{1!s}"),\n'.format(system.surface_volume_ratio.value, system.surface_volume_ratio.units))
+        elif isinstance(system, ConstantVIdealGasReactor) or isinstance(system, ConstantTPIdealGasReactor):
+            f.write('constantVIdealGasReactor(\n')
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    pressure = ' + format_pressure(system) + '\n')
+            f.write('    initialMoleFractions={\n')
+            for spcs, conc in system.initial_conditions.items():
+                if spcs in ['T', 'P']:
+                    continue
+                f.write('        "{0!s}": {1:g},\n'.format(spcs, conc))
+            f.write('    },\n')
+        elif isinstance(system, ConstantTVLiquidReactor):
+            f.write('constantTVLiquidReactor(\n')
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    initialConcentrations={\n')
+            for spcs, conc in system.initial_conditions.items():
+                if spcs in ['T', 'V']:
+                    continue
+                f.write('        "{0!s}": ({1:g},"{2!s}"),\n'.format(spcs, conc, 'mol/m^3'))
+            f.write('    },\n')
+        elif isinstance(system, MBSampledReactor):
+            f.write('mbsampledReactor(\n')
+            f.write('    temperature = ' + format_temperature(system) + '\n')
+            f.write('    pressure = ' + format_pressure(system) + '\n')
+            f.write('    initialMoleFractions={\n')
+            f.write(format_initial_mole_fractions(system))
+            f.write('    },\n')
+            f.write('    mbsamplingRate = ' + str(system.k_sampling.value_si) + ',\n')
+            f.write('    constantSpecies = ' + str([x.label for x in system.constantSpeciesList]) + ',\n')
         else:
             f.write('simpleReactor(\n')
             f.write('    temperature = ' + format_temperature(system) + '\n')
             f.write('    pressure = ' + format_pressure(system) + '\n')
             f.write('    initialMoleFractions={\n')
             f.write(format_initial_mole_fractions(system))
-        f.write('    },\n')
+            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):
@@ -1833,7 +1907,7 @@ def format_initial_mole_fractions(system):
     f.write('    maximumEdgeSpecies = {0:d},\n'.format(rmg.model_settings_list[0].maximum_edge_species))
     f.write('    minCoreSizeForPrune = {0:d},\n'.format(rmg.model_settings_list[0].min_core_size_for_prune))
     f.write('    minSpeciesExistIterationsForPrune = {0:d},\n'.format(rmg.model_settings_list[0].min_species_exist_iterations_for_prune))
-    f.write('    filterReactions = {0:d},\n'.format(rmg.model_settings_list[0].filter_reactions))
+    f.write('    filterReactions = {0},\n'.format(bool(rmg.model_settings_list[0].filter_reactions)))
     f.write('    filterThreshold = {0:g},\n'.format(rmg.model_settings_list[0].filter_threshold))
     f.write(')\n\n')
 
@@ -1903,7 +1977,7 @@ def formula(elements):
     f.write('    keepIrreversible = {0},\n'.format(rmg.keep_irreversible))
     f.write('    trimolecularProductReversible = {0},\n'.format(rmg.trimolecular_product_reversible))
     f.write('    verboseComments = {0},\n'.format(rmg.verbose_comments))
-    f.write('    wallTime = {0},\n'.format(rmg.walltime))
+    f.write('    wallTime = "{0}",\n'.format(rmg.walltime))
     f.write(')\n\n')
 
     f.close()
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index 3ce2ce4120..a2fe52feb1 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -807,8 +807,8 @@ def execute(self, initialize=True, **kwargs):
         self.done = False
 
         # determine min and max values for T and P (don't determine P values for liquid reactors)
-        self.Tmin = min([x.Trange[0].value_si if x.Trange else x.T.value_si for x in self.reaction_systems])
-        self.Tmax = max([x.Trange[1].value_si if x.Trange else x.T.value_si for x in self.reaction_systems])
+        self.Tmin = min([x.Trange[0].value_si if hasattr(x, "Trange") and x.Trange else x.T.value_si for x in self.reaction_systems])
+        self.Tmax = max([x.Trange[1].value_si if hasattr(x, "Trange") and x.Trange else x.T.value_si for x in self.reaction_systems])
         try:
             self.Pmin = min([x.Prange[0].value_si if hasattr(x, "Prange") and x.Prange else x.P.value_si for x in self.reaction_systems])
             self.Pmax = max([x.Prange[1].value_si if hasattr(x, "Prange") and x.Prange else x.P.value_si for x in self.reaction_systems])
diff --git a/rmgpy/solver/liquid.pyx b/rmgpy/solver/liquid.pyx
index 9ebc9889b6..fef047f215 100644
--- a/rmgpy/solver/liquid.pyx
+++ b/rmgpy/solver/liquid.pyx
@@ -91,8 +91,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
 
diff --git a/test/rmgpy/rmg/inputTest.py b/test/rmgpy/rmg/inputTest.py
index 00160d974b..6f5f3686f2 100644
--- a/test/rmgpy/rmg/inputTest.py
+++ b/test/rmgpy/rmg/inputTest.py
@@ -28,7 +28,7 @@
 ###############################################################################
 
 from unittest.mock import patch
-
+import rmgpy
 import rmgpy.rmg.input as inp
 from rmgpy.rmg.main import RMG
 from rmgpy.rmg.model import CoreEdgeReactionModel
@@ -447,3 +447,538 @@ def test_completed_networks_none(self):
         
         # Check that no networks were added
         assert len(rmg.reaction_model.completed_pdep_networks) == 0
+
+
+class TestWriteInputFile:
+    """
+    Contains unit test for writing input files for each of the reactor types:
+
+        'simpleReactor': simple_reactor, ✅
+        'constantVIdealGasReactor' : constant_V_ideal_gas_reactor, ✅
+        'constantTPIdealGasReactor' : constant_TP_ideal_gas_reactor, ✅
+        'liquidSurfaceReactor' : liquid_cat_reactor, ✅
+        'constantTVLiquidReactor': constant_T_V_liquid_reactor,
+        'liquidReactor': liquid_reactor, ✅
+        'surfaceReactor': surface_reactor, ✅
+        'mbsampledReactor': mb_sampled_reactor,
+
+    """
+    def setup_method(self):
+        """This method is run before every test in this class"""
+        global rmg
+        rmg.reaction_systems = []
+
+    def test_write_superminimal_input(self):
+        """
+        Test that we can write superminimal input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        superminimal_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/superminimal/input.py')
+        superminimal_output_file = 'temp_superminimal_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(superminimal_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P.value_si
+        initialMoleFractions = {k.label: v for k, v in rmg.reaction_systems[0].initial_mole_fractions.items()}
+        for term in rmg.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                termination_time = term.time.value_si
+            elif hasattr(term, 'conversion'):
+                termination_conversion = term.conversion
+                termination_converstion_species = term.species.label
+
+        inp.save_input_file(superminimal_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(superminimal_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P.value_si == P
+        output_mol_fractions = {k.label: v for k, v in rmg1.reaction_systems[0].initial_mole_fractions.items()}
+        assert output_mol_fractions == initialMoleFractions
+        for term in rmg1.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                assert term.time.value_si == termination_time
+            elif hasattr(term, 'conversion'):
+                assert term.conversion == termination_conversion
+                assert term.species.label == termination_converstion_species
+
+        # clean up
+        os.remove(superminimal_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_superminimal_and_run(self):
+        """
+        Test that we can write superminimal input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        superminimal_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/superminimal/input.py')
+        new_run_dir = 'temp_superminimal_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        superminimal_output_file = os.path.join(new_run_dir, 'temp_superminimal_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(superminimal_input_file, rmg)
+        inp.save_input_file(superminimal_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), superminimal_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    def test_write_min_surf_input(self):
+        """
+        Test that we can write the minimal surface input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        min_surf_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/minimal_surface/input.py')
+        min_surf_output_file = 'temp_min_surf_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(min_surf_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P_initial.value_si
+        initialMoleFractions = {k.label: v for k, v in rmg.reaction_systems[0].initial_gas_mole_fractions.items()}
+        initialSurfaceCoverages = {k.label: v for k, v in rmg.reaction_systems[0].initial_surface_coverages.items()}
+        for term in rmg.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                termination_time = term.time.value_si
+            elif hasattr(term, 'conversion'):
+                termination_conversion = term.conversion
+                termination_converstion_species = term.species.label
+            elif hasattr(term, 'ratio'):
+                termination_ratio = term.ratio
+
+        binding_energies = {k: v.value_si for k, v in rmg.binding_energies.items()}
+        surface_site_density = rmg.surface_site_density.value_si
+
+        inp.save_input_file(min_surf_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(min_surf_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P_initial.value_si == P
+        output_mol_fractions = {k.label: v for k, v in rmg1.reaction_systems[0].initial_gas_mole_fractions.items()}
+        assert output_mol_fractions == initialMoleFractions
+        output_surface_coverages = {k.label: v for k, v in rmg1.reaction_systems[0].initial_surface_coverages.items()}
+        assert output_surface_coverages == initialSurfaceCoverages
+
+        output_binding_energies = {k: v.value_si for k, v in rmg1.binding_energies.items()}
+        assert output_binding_energies == binding_energies
+
+        assert rmg1.surface_site_density.value_si == surface_site_density
+
+        for term in rmg1.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                assert term.time.value_si == termination_time
+            elif hasattr(term, 'conversion'):
+                assert term.conversion == termination_conversion
+                assert term.species.label == termination_converstion_species
+            elif hasattr(term, 'ratio'):
+                assert term.ratio == termination_ratio
+
+        # clean up
+        os.remove(min_surf_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_min_surf_and_run(self):
+        """
+        Test that we can write minimal surface input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        min_surf_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/minimal_surface/input.py')
+        new_run_dir = 'temp_min_surf_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        min_surf_output_file = os.path.join(new_run_dir, 'temp_min_surf_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(min_surf_input_file, rmg)
+        inp.save_input_file(min_surf_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), min_surf_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    @pytest.mark.skip(reason="Slow because it has to compile Julia")
+    def test_write_liquid_cat_input(self):
+        """
+        Test that we can write liquid catalyst input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        liquid_cat_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/liquid_cat/input.py')
+        liquid_cat_output_file = 'temp_liquid_cat_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(liquid_cat_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        tf = rmg.reaction_systems[0].tf
+        liquid_init_conditions = {k: v for k, v in rmg.reaction_systems[0].initial_conditions['liquid'].items()}
+        surf_init_conditions = {k: v for k, v in rmg.reaction_systems[0].initial_conditions['surface'].items()}
+
+        termination_species = rmg.reaction_systems[0].terminations[0][0].label
+        termination_conversion = rmg.reaction_systems[0].terminations[0][1]
+
+        inp.save_input_file(liquid_cat_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(liquid_cat_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].tf == tf
+
+        liquid_init_conditions_output = {k: v for k, v in rmg1.reaction_systems[0].initial_conditions['liquid'].items()}
+        assert pytest.approx(liquid_init_conditions_output) == liquid_init_conditions
+        surf_init_conditions_output = {k: v for k, v in rmg1.reaction_systems[0].initial_conditions['surface'].items()}
+        assert pytest.approx(surf_init_conditions_output) == surf_init_conditions
+
+        assert rmg1.reaction_systems[0].terminations[0][0].label == termination_species
+        assert rmg1.reaction_systems[0].terminations[0][1] == termination_conversion
+
+        # clean up
+        os.remove(liquid_cat_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_liquid_cat_and_run(self):
+        """
+        Test that we can write liquid catalyst input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        liquid_cat_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/liquid_cat/input.py')
+        new_run_dir = 'temp_liquid_cat_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        liquid_cat_output_file = os.path.join(new_run_dir, 'temp_liquid_cat_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(liquid_cat_input_file, rmg)
+        inp.save_input_file(liquid_cat_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), '-t', '00:00:01:30', liquid_cat_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    def test_write_liquid_input(self):
+        """
+        Test that we can write the liquid reactor input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        liquid_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/liquid_phase/input.py')
+        liquid_output_file = 'temp_liquid_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(liquid_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P.value_si
+        initialConcentrations = {k.label: v for k, v in rmg.reaction_systems[0].initial_concentrations.items()}
+        for term in rmg.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                termination_time = term.time.value_si
+        solvent = rmg.solvent
+
+        inp.save_input_file(liquid_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(liquid_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P.value_si == P
+        output_concentrations = {k.label: v for k, v in rmg1.reaction_systems[0].initial_concentrations.items()}
+        assert pytest.approx(output_concentrations) == initialConcentrations
+        for term in rmg1.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                assert term.time.value_si == termination_time
+        assert rmg1.solvent == solvent
+
+        # clean up
+        os.remove(liquid_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_liquid_and_run(self):
+        """
+        Test that we can write liquid reactor input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        liquid_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/liquid_phase/input.py')
+        new_run_dir = 'temp_liquid_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        liquid_output_file = os.path.join(new_run_dir, 'temp_liquid_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(liquid_input_file, rmg)
+        inp.save_input_file(liquid_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), '-t', '00:00:01:30', liquid_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    @pytest.mark.skip(reason="Slow because it has to compile Julia")
+    def test_write_constantVIdealGasReactor(self):
+        """
+        Test that we can write constant volume ideal gas reactor input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        rms_constant_V_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/rms_constant_V/input.py')
+        rms_constant_V_output_file = 'temp_rms_constant_V_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(rms_constant_V_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P.value_si
+        tf = rmg.reaction_systems[0].tf
+        init_conditions = {k: v for k, v in rmg.reaction_systems[0].initial_conditions.items()}
+
+        termination_species = rmg.reaction_systems[0].terminations[0][0].label
+        termination_conversion = rmg.reaction_systems[0].terminations[0][1]
+        termination_time = rmg.reaction_systems[0].terminations[1].time
+
+        inp.save_input_file(rms_constant_V_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(rms_constant_V_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P.value_si == P
+        assert rmg1.reaction_systems[0].tf == tf
+
+        new_init_conditions = {k: v for k, v in rmg1.reaction_systems[0].initial_conditions.items()}
+        assert pytest.approx(new_init_conditions) == init_conditions
+
+        assert rmg1.reaction_systems[0].terminations[0][0].label == termination_species
+        assert rmg1.reaction_systems[0].terminations[0][1] == termination_conversion
+        assert rmg1.reaction_systems[0].terminations[1].time == termination_time
+
+        # clean up
+        os.remove(rms_constant_V_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_constantVIdealGasReactor_and_run(self):
+        """
+        Test that we can write constant volume ideal gas reactor input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        constant_V_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/rms_constant_V/input.py')
+        new_run_dir = 'temp_constant_V_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        constant_V_output_file = os.path.join(new_run_dir, 'temp_constant_V_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(constant_V_input_file, rmg)
+        inp.save_input_file(constant_V_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), '-t', '00:00:01:30', constant_V_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    @pytest.mark.skip(reason="Slow because it has to compile Julia")
+    def test_write_constantTPdealGasReactor(self):
+        """
+        Test that we can write constant TP ideal gas reactor input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        rms_constant_TP_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/nox_transitory_edge/input.py')
+        rms_constant_TP_output_file = 'temp_constant_TP_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(rms_constant_TP_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P.value_si
+        tf = rmg.reaction_systems[0].tf
+        init_conditions = {k: v for k, v in rmg.reaction_systems[0].initial_conditions.items()}
+
+        termination_species = rmg.reaction_systems[0].terminations[0][0].label
+        termination_conversion = rmg.reaction_systems[0].terminations[0][1]
+        termination_time = rmg.reaction_systems[0].terminations[1].time
+
+        inp.save_input_file(rms_constant_TP_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(rms_constant_TP_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P.value_si == P
+        assert rmg1.reaction_systems[0].tf == tf
+
+        new_init_conditions = {k: v for k, v in rmg1.reaction_systems[0].initial_conditions.items()}
+        assert pytest.approx(new_init_conditions, rel=1e-4) == init_conditions
+
+        assert rmg1.reaction_systems[0].terminations[0][0].label == termination_species
+        assert rmg1.reaction_systems[0].terminations[0][1] == termination_conversion
+        assert rmg1.reaction_systems[0].terminations[1].time == termination_time
+
+        # clean up
+        os.remove(rms_constant_TP_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_constantTPIdealGasReactor_and_run(self):
+        """
+        Test that we can write constant TP ideal gas reactor input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        constant_TP_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../examples/rmg/nox_transitory_edge/input.py')
+        new_run_dir = 'temp_constant_TP_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        constant_TP_output_file = os.path.join(new_run_dir, 'temp_constant_TP_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(constant_TP_input_file, rmg)
+        inp.save_input_file(constant_TP_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), '-t', '00:00:01:30', constant_TP_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    @pytest.mark.skip(reason="Slow because it has to compile Julia")
+    def test_write_constantTVLiquidReactor(self):
+        """
+        Test that we can write constant TV liquid reactor input file and read it back in with the same values.
+        """
+        import os
+        import rmgpy
+        rms_constant_TV_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../test/regression/RMS_CSTR_liquid_oxidation/input.py')
+        rms_constant_TV_output_file = 'temp_constant_TV_liquid_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(rms_constant_TV_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        V = rmg.reaction_systems[0].initial_conditions['V']
+        tf = rmg.reaction_systems[0].tf
+        init_conditions = {k: v for k, v in rmg.reaction_systems[0].initial_conditions.items()}
+
+        termination_species = rmg.reaction_systems[0].terminations[0][0].label
+        termination_conversion = rmg.reaction_systems[0].terminations[0][1]
+        termination_time = rmg.reaction_systems[0].terminations[1].time
+
+        inp.save_input_file(rms_constant_TV_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(rms_constant_TV_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].initial_conditions['V'] == V
+        assert rmg1.reaction_systems[0].tf == tf
+
+        new_init_conditions = {k: v for k, v in rmg1.reaction_systems[0].initial_conditions.items()}
+        assert pytest.approx(new_init_conditions, rel=1e-4) == init_conditions
+
+        assert rmg1.reaction_systems[0].terminations[0][0].label == termination_species
+        assert rmg1.reaction_systems[0].terminations[0][1] == termination_conversion
+        assert rmg1.reaction_systems[0].terminations[1].time == termination_time
+
+        # clean up
+        os.remove(rms_constant_TV_output_file)
+
+    @pytest.mark.skip(reason="Slow test that runs a full RMG job")
+    def test_write_constantTVLiquidReactor_and_run(self):
+        """
+        Test that we can write constant TV liquid reactor input file and then run RMG without errors
+        """
+        import os
+        import rmgpy
+        import shutil
+
+        constant_TV_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../test/regression/RMS_CSTR_liquid_oxidation/input.py')
+        new_run_dir = 'temp_constant_TV_run'
+        os.makedirs(new_run_dir, exist_ok=True)
+        constant_TV_output_file = os.path.join(new_run_dir, 'temp_constant_TV_input.py')
+
+        rmg = RMG()
+        inp.read_input_file(constant_TV_input_file, rmg)
+        inp.save_input_file(constant_TV_output_file, rmg)
+
+        # run RMG with the new input file
+        import subprocess
+        subprocess.run(['python', os.path.join(rmgpy.settings['test_data.directory'], '../../../rmg.py'), '-t', '00:00:01:30', constant_TV_output_file], check=True)
+
+        # clean up
+        shutil.rmtree(new_run_dir)
+
+    def test_MBSampledReactor_write(self):
+        """
+        Test that we can write MB sampled reactor input file and read it back in with the same values.
+        Note that the MBSampledReactor is not intended to be used with a standard RMG job, so there's no point in running it as a test
+        """
+        import os
+        import rmgpy
+        mbsampled_input_file = os.path.join(rmgpy.settings['test_data.directory'], '../../../rmgpy/tools/data/sim/mbSampled/input.py')
+        mbsampled_output_file = 'temp_mbsampled_input.py'
+
+        rmg = RMG()
+        inp.read_input_file(mbsampled_input_file, rmg)
+
+        # read a bunch of values in from input file to check they are the same after writing
+        T = rmg.reaction_systems[0].T.value_si
+        P = rmg.reaction_systems[0].P.value_si
+        sampling_rate = rmg.reaction_systems[0].k_sampling.value_si
+
+        initialMoleFractions = {k.label: v for k, v in rmg.reaction_systems[0].initial_mole_fractions.items()}
+
+        for term in rmg.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                termination_time = term.time.value_si
+
+        inp.save_input_file(mbsampled_output_file, rmg)
+        # read it back in and confirm all the values match
+        rmg1 = RMG()
+        inp.read_input_file(mbsampled_output_file, rmg1)
+        assert rmg1.reaction_systems[0].T.value_si == T
+        assert rmg1.reaction_systems[0].P.value_si == P
+        assert rmg1.reaction_systems[0].k_sampling.value_si == sampling_rate
+
+        new_initialMoleFractions = {k.label: v for k, v in rmg1.reaction_systems[0].initial_mole_fractions.items()}
+        assert pytest.approx(new_initialMoleFractions, rel=1e-4) == initialMoleFractions
+
+        for term in rmg1.reaction_systems[0].termination:
+            if hasattr(term, 'time'):
+                assert term.time.value_si == termination_time
+
+        # clean up
+        os.remove(mbsampled_output_file)