diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 7b5a39c162..44ef22a240 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -122,6 +122,12 @@ jobs:
           activate-environment: rmg_env
           use-mamba: true
 
+      # installs the extra RMS conda dependencies
+      - name: Add RMS dependencies
+        run: |
+          mamba install -c conda-forge julia=1.9.1 pyjulia>=0.6 
+          mamba install -c rmg pyrms diffeqpy
+
       # list the environment for debugging purposes
       - name: mamba info
         run: |
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 8ef810f67f..1da6623fdc 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -36,6 +36,12 @@ jobs:
           activate-environment: rmg_env
           use-mamba: true
 
+      # installs the extra RMS conda dependencies
+      - name: Add RMS dependencies
+        run: |
+          mamba install -c conda-forge julia=1.9.1 pyjulia>=0.6 
+          mamba install -c rmg pyrms diffeqpy
+
       - name: Install sphinx
         run: mamba install -y sphinx
 
diff --git a/Dockerfile b/Dockerfile
index a263199c89..d34adc2a2c 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -28,10 +28,6 @@ RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh &
     rm Miniconda3-latest-Linux-x86_64.sh
 ENV PATH="$PATH:/miniconda/bin"
 
-# Set solver backend to mamba for speed
-RUN conda install -n base conda-libmamba-solver && \
-    conda config --set solver libmamba
-
 # Set Bash as the default shell for following commands
 SHELL ["/bin/bash", "-c"]
 
@@ -44,8 +40,7 @@ RUN git clone --single-branch --branch main --depth 1 https://github.com/Reactio
 
 WORKDIR /rmg/RMG-Py
 # build the conda environment
-RUN conda env create --file environment.yml && \
-    conda clean --all --yes
+RUN conda env create --file environment.yml
 
 # This runs all subsequent commands inside the rmg_env conda environment
 #
@@ -54,6 +49,10 @@ RUN conda env create --file environment.yml && \
 # in a Dockerfile build script)
 SHELL ["conda", "run", "--no-capture-output", "-n", "rmg_env", "/bin/bash", "-c"]
 
+RUN conda install -c conda-forge julia=1.9.1 pyjulia>=0.6 && \
+    conda install -c rmg pyrms diffeqpy && \
+    conda clean --all --yes
+
 # Set environment variables as directed in the RMG installation instructions
 ENV RUNNER_CWD=/rmg
 ENV PYTHONPATH="$RUNNER_CWD/RMG-Py:$PYTHONPATH"
diff --git a/documentation/source/users/rmg/installation/anacondaDeveloper.rst b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
index 4a79244d0c..b081bab115 100644
--- a/documentation/source/users/rmg/installation/anacondaDeveloper.rst
+++ b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
@@ -25,6 +25,11 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
 
       Note that you should reinitialize or restart your terminal in order for the changes to take effect, as the installer will tell you.
 
+#. If your `conda` version is older than 23.10.0, switch the solver backend to `libmamba` ::
+
+    conda install -n base conda-libmamba-solver
+    conda config --set solver libmamba
+
 #. There are a few system-level dependencies which are required and should not be installed via Conda. These include
    `Git <https://git-scm.com/>`_ for version control, `GNU Make <https://www.gnu.org/software/make/>`_, and the C and C++ compilers from the `GNU Compiler Collection (GCC) <https://gcc.gnu.org/>`_ for compiling RMG.
 
@@ -71,11 +76,6 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
 
    For information on using ``ssh`` with GitHub see the `Connecting to GitHub with SSH <https://docs.github.com/en/authentication/connecting-to-github-with-ssh>`_
 
-#. Switch the conda solver backend to speed up creation of the RMG environment ::
-
-    conda install -n base conda-libmamba-solver
-    conda config --set solver libmamba
-
 #. Navigate to the RMG-Py directory ::
 
     cd RMG-Py
@@ -110,16 +110,11 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
 
     conda activate rmg_env
 
-#. Switch the conda solver to libmamba again, to accelerate any changes you might make to this conda environment in the future::
-
-    conda config --set solver libmamba
-
 #. Compile RMG-Py after activating the conda environment ::
 
     make
 
 #. Modify environment variables. Add RMG-Py to the PYTHONPATH to ensure that you can access RMG modules from any folder.
-   *This is important before the next step in which julia dependencies are installed.*
    Also, add your RMG-Py folder to PATH to launch ``rmg.py`` from any folder.
 
    In general, these commands should be placed in the appropriate shell initialization file.
@@ -134,16 +129,21 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
 
    Be sure to either close and reopen your terminal to refresh your environment variables (``source ~/.bashrc`` or ``source ~/.zshrc``).
 
-#. Install and Link Julia dependencies: ::
+#. **Optional (Recommended)**: Install and Link Julia dependencies. Ensure that you have modified your environment variables as described above, and then run the following: ::
 
      julia -e 'using Pkg; Pkg.add("PyCall");Pkg.build("PyCall");Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator;'
 
      python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
 
+    Installing these dependencies will allow using ``method='ode'`` when solving the Master Equation with Arkane and using ``ReactionMechanismSimulator.jl``-based reactors in RMG.
 
 #. Finally, you can run RMG from any location by typing the following (given that you have prepared the input file as ``input.py`` in the current folder). ::
 
-    python-jl replace/with/path/to/rmg.py input.py
+     python-jl rmg.py input.py
+    
+    or, if the Julia dependencies are not installed: ::
+
+     python rmg.py input.py
 
 You may now use RMG-Py, Arkane, as well as any of the :ref:`Standalone Modules <modules>` included in the RMG-Py package.
 
diff --git a/environment.yml b/environment.yml
index 48a4da34a3..ac5bd4b6e4 100644
--- a/environment.yml
+++ b/environment.yml
@@ -16,6 +16,7 @@
 #   made dependency list more explicit (@JacksonBurns).
 # - October 16, 2023 Switched RDKit and descripatastorus to conda-forge,
 #   moved diffeqpy to pip and (temporarily) removed chemprop
+# - Mar 11, 2024 Removed Julia dependencies, now considered optional
 #
 name: rmg_env
 channels:
@@ -48,10 +49,6 @@ dependencies:
   - conda-forge::openbabel >= 3
   - conda-forge::rdkit >=2022.09.1
 
-# general-purpose external software tools
-  - conda-forge::julia=1.9.1
-  - conda-forge::pyjulia >=0.6
-
 # Python tools
   - python >=3.7
   - coverage
@@ -88,21 +85,24 @@ dependencies:
 # packages we maintain
   - rmg::pydas >=1.0.3
   - rmg::pydqed >=1.0.3
-  - rmg::pyrms
   - rmg::symmetry
 
-# packages we would like to stop maintaining (and why)
-  - rmg::diffeqpy
-  # we should use the official verison https://github.com/SciML/diffeqpy),
-  # rather than ours (which is only made so that we can get it from conda)
-  # It is only on pip, so we will need to do something like:
-  # https://stackoverflow.com/a/35245610
-  # Note that _some other_ dep. in this list requires diffeqpy in its recipe
-  # which will cause it to be downloaded from the rmg conda channel
-
 # conda mutex metapackage
   - nomkl
 
+# optional dependencies for using ReactionMechanismSimulator
+# remove the leading '#' to install the required dependencies
+  # - conda-forge::julia=1.9.1
+  # - conda-forge::pyjulia >=0.6
+  # - rmg::pyrms
+  # - rmg::diffeqpy
+
+# Note about diffeqpy:
+# we should use the official verison https://github.com/SciML/diffeqpy),
+# rather than ours (which is only made so that we can get it from conda)
+# It is only on pip, so we will need to do something like:
+# https://stackoverflow.com/a/35245610
+
 # additional packages that are required, but not specified here (and why)
   # pydqed, pydas, mopac, and likely others require a fortran compiler (specifically gfortran)
   # in the environment. Normally we would add this to the environment file with
diff --git a/rmgpy/pdep/sls.py b/rmgpy/pdep/sls.py
index b183a0bba8..a292e8391a 100644
--- a/rmgpy/pdep/sls.py
+++ b/rmgpy/pdep/sls.py
@@ -32,9 +32,8 @@
 and implementing the SLS master equation reduction method
 """
 import sys
+import logging
 
-from diffeqpy import de
-from julia import Main
 import scipy.sparse as sparse
 import numpy as np
 import scipy.linalg
@@ -46,6 +45,17 @@
 from rmgpy.pdep.me import generate_full_me_matrix, states_to_configurations
 from rmgpy.statmech.translation import IdealGasTranslation
 
+NO_JULIA = False
+try:
+    from diffeqpy import de
+    from julia import Main
+except Exception as e:
+    logging.info(
+        f"Unable to import Julia dependencies, original error: {str(e)}"
+        ". Master equation method 'ode' will not be available on this execution."
+    )
+    NO_JULIA = True
+
 
 def get_initial_condition(network, x0, indices):
     """
@@ -150,6 +160,11 @@ def get_rate_coefficients_SLS(network, T, P, method="mexp", neglect_high_energy_
     tau = np.abs(1.0 / fastest_reaction)
 
     if method == "ode":
+        if NO_JULIA:
+            raise RuntimeError(
+                "Required Julia dependencies for method 'ode' are not installed.\n"
+                "Please follow the steps to install Julia dependencies at https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html."
+            )
         f = Main.eval(
             """
 function f(u,M,t)
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index e04cd82e2e..ba9b86f87a 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -48,9 +48,10 @@
 from rmgpy.solver.surface import SurfaceReactor
 from rmgpy.util import as_list
 from rmgpy.data.surface import MetalDatabase
-from rmgpy.rmg.reactors import Reactor, ConstantVIdealGasReactor, ConstantTLiquidSurfaceReactor, ConstantTVLiquidReactor, ConstantTPIdealGasReactor
 from rmgpy.data.vaporLiquidMassTransfer import liquidVolumetricMassTransferCoefficientPowerLaw
 from rmgpy.molecule.fragment import Fragment
+from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor, ConstantVIdealGasReactor, ConstantTLiquidSurfaceReactor, ConstantTVLiquidReactor, ConstantTPIdealGasReactor, NO_JULIA
+ 
 
 ################################################################################
 
@@ -1558,6 +1559,11 @@ def read_input_file(path, rmg0):
         exec(f.read(), global_context, local_context)
     except (NameError, TypeError, SyntaxError) as e:
         logging.error('The input file "{0}" was invalid:'.format(full_path))
+        if NO_JULIA:
+            logging.error(
+                "During runtime, import of Julia dependencies failed. To use phase systems and RMS reactors, install RMG-Py with RMS."
+                " (https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html)"
+            )
         logging.exception(e)
         raise
     finally:
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index 8df771980d..8a9ebf2d2a 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -80,7 +80,7 @@
 from rmgpy.tools.plot import plot_sensitivity
 from rmgpy.tools.uncertainty import Uncertainty, process_local_results
 from rmgpy.yml import RMSWriter
-from rmgpy.rmg.reactors import Reactor
+from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor
 
 ################################################################################
 
@@ -507,6 +507,12 @@ def initialize(self, **kwargs):
 
         # Read input file
         self.load_input(self.input_file)
+        
+        # Check if ReactionMechanismSimulator reactors are being used
+        # if RMS is not installed but the user attempted to use it, the load_input_file would have failed
+        # if RMS is installed but they did not use it, we can avoid extra work
+        # if RMS is not installed and they did not use it, we avoid calling certain functions that would raise an error
+        requires_rms = any(isinstance(reactor_system, RMSReactor) for reactor_system in self.reaction_systems)
 
         if kwargs.get("restart", ""):
             import rmgpy.rmg.input
@@ -550,10 +556,10 @@ def initialize(self, **kwargs):
         self.load_database()
 
         for spec in self.initial_species:
-            self.reaction_model.add_species_to_edge(spec)
+            self.reaction_model.add_species_to_edge(spec, requires_rms=requires_rms)
 
         for reaction_system in self.reaction_systems:
-            if isinstance(reaction_system, Reactor):
+            if isinstance(reaction_system, RMSReactor):
                 reaction_system.finish_termination_criteria()
 
         # Load restart seed mechanism (if specified)
@@ -618,12 +624,12 @@ def initialize(self, **kwargs):
         # Seed mechanisms: add species and reactions from seed mechanism
         # DON'T generate any more reactions for the seed species at this time
         for seed_mechanism in self.seed_mechanisms:
-            self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False)
+            self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False, requires_rms=requires_rms)
 
         # Reaction libraries: add species and reactions from reaction library to the edge so
         # that RMG can find them if their rates are large enough
         for library, option in self.reaction_libraries:
-            self.reaction_model.add_reaction_library_to_edge(library)
+            self.reaction_model.add_reaction_library_to_edge(library, requires_rms=requires_rms)
 
         # Also always add in a few bath gases (since RMG-Java does)
         for label, smiles in [("Ar", "[Ar]"), ("He", "[He]"), ("Ne", "[Ne]"), ("N2", "N#N")]:
@@ -695,35 +701,37 @@ def initialize(self, **kwargs):
         # This is necessary so that the PDep algorithm can identify the bath gas
         for spec in self.initial_species:
             if not spec.reactive:
-                self.reaction_model.enlarge(spec)
+                self.reaction_model.enlarge(spec, requires_rms=requires_rms)
         for spec in self.initial_species:
             if spec.reactive:
-                self.reaction_model.enlarge(spec)
+                self.reaction_model.enlarge(spec, requires_rms=requires_rms)
 
         # chatelak: store constant SPC indices in the reactor attributes if any constant SPC provided in the input file
         # advantages to write it here: this is run only once (as species indexes does not change over the generation)
         if self.solvent is not None:
             for index, reaction_system in enumerate(self.reaction_systems):
                 if (
-                    not isinstance(reaction_system, Reactor) and reaction_system.const_spc_names is not None
-                ):  # if no constant species provided do nothing
+                    not isinstance(reaction_system, RMSReactor)
+                ) and reaction_system.const_spc_names is not None:  # if no constant species provided do nothing
                     reaction_system.get_const_spc_indices(self.reaction_model.core.species)  # call the function to identify indices in the solver
 
         self.initialize_reaction_threshold_and_react_flags()
         if self.filter_reactions and self.init_react_tuples:
-            self.react_init_tuples()
+            self.react_init_tuples(requires_rms=requires_rms)
         self.reaction_model.initialize_index_species_dict()
 
         self.initialize_seed_mech()
+        return requires_rms
 
-    def register_listeners(self):
+    def register_listeners(self, requires_rms=False):
         """
         Attaches listener classes depending on the options
         found in the RMG input file.
         """
 
         self.attach(ChemkinWriter(self.output_directory))
-        self.attach(RMSWriter(self.output_directory))
+        if requires_rms:
+            self.attach(RMSWriter(self.output_directory))
 
         if self.generate_output_html:
             self.attach(OutputHTMLWriter(self.output_directory))
@@ -735,7 +743,7 @@ def register_listeners(self):
 
         if self.save_simulation_profiles:
             for index, reaction_system in enumerate(self.reaction_systems):
-                if isinstance(reaction_system, Reactor):
+                if requires_rms and isinstance(reaction_system, RMSReactor):
                     typ = type(reaction_system)
                     raise InputError(f"save_simulation_profiles=True not compatible with reactor of type {typ}")
                 reaction_system.attach(SimulationProfileWriter(self.output_directory, index, self.reaction_model.core.species))
@@ -749,10 +757,10 @@ def execute(self, initialize=True, **kwargs):
         """
 
         if initialize:
-            self.initialize(**kwargs)
+            requires_rms = self.initialize(**kwargs)
 
         # register listeners
-        self.register_listeners()
+        self.register_listeners(requires_rms=requires_rms)
 
         self.done = False
 
@@ -779,7 +787,7 @@ def execute(self, initialize=True, **kwargs):
             # Update react flags
             if self.filter_reactions:
                 # Run the reaction system to update threshold and react flags
-                if isinstance(reaction_system, Reactor):
+                if requires_rms and isinstance(reaction_system, RMSReactor):
                     self.update_reaction_threshold_and_react_flags(
                         rxn_sys_unimol_threshold=np.zeros((len(self.reaction_model.core.species),), bool),
                         rxn_sys_bimol_threshold=np.zeros((len(self.reaction_model.core.species), len(self.reaction_model.core.species)), bool),
@@ -822,6 +830,7 @@ def execute(self, initialize=True, **kwargs):
                 unimolecular_react=self.unimolecular_react,
                 bimolecular_react=self.bimolecular_react,
                 trimolecular_react=self.trimolecular_react,
+                requires_rms=requires_rms,
             )
 
         if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
@@ -834,7 +843,7 @@ def execute(self, initialize=True, **kwargs):
             )
 
         if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
-            self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species)
+            self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species, requires_rms=requires_rms)
 
         logging.info("Completed initial enlarge edge step.\n")
 
@@ -900,7 +909,7 @@ def execute(self, initialize=True, **kwargs):
                             prune = False
 
                         try:
-                            if isinstance(reaction_system, Reactor):
+                            if requires_rms and isinstance(reaction_system, RMSReactor):
                                 (
                                     terminated,
                                     resurrected,
@@ -993,7 +1002,7 @@ def execute(self, initialize=True, **kwargs):
 
                         # Add objects to enlarge to the core first
                         for objectToEnlarge in objects_to_enlarge:
-                            self.reaction_model.enlarge(objectToEnlarge)
+                            self.reaction_model.enlarge(objectToEnlarge, requires_rms=requires_rms)
 
                         if model_settings.filter_reactions:
                             # Run a raw simulation to get updated reaction system threshold values
@@ -1002,7 +1011,7 @@ def execute(self, initialize=True, **kwargs):
                             temp_model_settings.tol_keep_in_edge = 0
                             if not resurrected:
                                 try:
-                                    if isinstance(reaction_system, Reactor):
+                                    if requires_rms and isinstance(reaction_system, RMSReactor):
                                         (
                                             terminated,
                                             resurrected,
@@ -1071,7 +1080,7 @@ def execute(self, initialize=True, **kwargs):
                                     skip_update=True,
                                 )
                                 logging.warning(
-                                    "Reaction thresholds/flags for Reaction System {0} was not updated due " "to resurrection".format(index + 1)
+                                    "Reaction thresholds/flags for Reaction System {0} was not updated due to resurrection".format(index + 1)
                                 )
 
                             logging.info("")
@@ -1094,13 +1103,14 @@ def execute(self, initialize=True, **kwargs):
                             unimolecular_react=self.unimolecular_react,
                             bimolecular_react=self.bimolecular_react,
                             trimolecular_react=self.trimolecular_react,
+                            requires_rms=requires_rms,
                         )
 
                         if old_edge_size != len(self.reaction_model.edge.reactions) or old_core_size != len(self.reaction_model.core.reactions):
                             reactor_done = False
 
                         if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
-                            self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species)
+                            self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species, requires_rms=requires_rms)
 
                         max_num_spcs_hit = len(self.reaction_model.core.species) >= model_settings.max_num_species
 
@@ -1127,6 +1137,7 @@ def execute(self, initialize=True, **kwargs):
                             model_settings.tol_move_to_core,
                             model_settings.maximum_edge_species,
                             model_settings.min_species_exist_iterations_for_prune,
+                            requires_rms=requires_rms,
                         )
                         # Perform garbage collection after pruning
                         collected = gc.collect()
@@ -1891,7 +1902,7 @@ def initialize_reaction_threshold_and_react_flags(self):
                         if self.trimolecular:
                             self.trimolecular_react[:num_restart_spcs, :num_restart_spcs, :num_restart_spcs] = False
 
-    def react_init_tuples(self):
+    def react_init_tuples(self, requires_rms=False):
         """
         Reacts tuples given in the react block
         """
@@ -1924,6 +1935,7 @@ def react_init_tuples(self):
             unimolecular_react=self.unimolecular_react,
             bimolecular_react=self.bimolecular_react,
             trimolecular_react=self.trimolecular_react,
+            requires_rms=requires_rms,
         )
 
     def update_reaction_threshold_and_react_flags(
@@ -2218,7 +2230,7 @@ def __init__(self, reaction_system, bspc):
                 if isinstance(value, list):
                     self.Ranges[key] = [v.value_si for v in value]
 
-        if isinstance(reaction_system, Reactor):
+        if isinstance(reaction_system, RMSReactor):
             self.tmax = reaction_system.tf
         else:
             for term in reaction_system.termination:
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 1ab45dc67c..ca67effe89 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -57,7 +57,8 @@
 from rmgpy.species import Species
 from rmgpy.thermo.thermoengine import submit
 from rmgpy.rmg.decay import decay_species
-from rmgpy.rmg.reactors import PhaseSystem, Phase, Interface, Reactor
+from rmgpy.rmg.reactionmechanismsimulator_reactors import PhaseSystem, Phase, Interface, NO_JULIA
+from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor
 from rmgpy.molecule.fragment import Fragment
 
 ################################################################################
@@ -75,7 +76,7 @@ def __init__(self, species=None, reactions=None, phases=None, interfaces={}):
         if phases is None:
             phases = {"Default": Phase(), "Surface": Phase()}
             interfaces = {frozenset({"Default", "Surface"}): Interface(list(phases.values()))}
-        self.phase_system = PhaseSystem(phases, interfaces)
+        self.phase_system = None if NO_JULIA else PhaseSystem(phases, interfaces)
 
     def __reduce__(self):
         """
@@ -349,9 +350,10 @@ def make_new_species(self, object, label="", reactive=True, check_existing=True,
         orilabel = spec.label
         label = orilabel
         i = 2
-        while any([label in phase.names for phase in self.edge.phase_system.phases.values()]):
-            label = orilabel + "-" + str(i)
-            i += 1
+        if self.edge.phase_system:  # None when RMS not installed
+            while any([label in phase.names for phase in self.edge.phase_system.phases.values()]):
+                label = orilabel + "-" + str(i)
+                i += 1
         spec.label = label
 
         logging.debug("Creating new species %s", spec.label)
@@ -589,7 +591,7 @@ def make_new_pdep_reaction(self, forward):
 
         return forward
 
-    def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None):
+    def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None, requires_rms=False):
         """
         Enlarge a reaction model by processing the objects in the list `new_object`.
         If `new_object` is a
@@ -636,7 +638,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
 
                 # Add new species
                 if new_species not in self.core.species:
-                    reactions_moved_from_edge = self.add_species_to_core(new_species)
+                    reactions_moved_from_edge = self.add_species_to_core(new_species, requires_rms=requires_rms)
                 else:
                     reactions_moved_from_edge = []
 
@@ -644,7 +646,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
                 pdep_network, new_species = new_object
                 new_reactions.extend(pdep_network.explore_isomer(new_species))
 
-                self.process_new_reactions(new_reactions, new_species, pdep_network)
+                self.process_new_reactions(new_reactions, new_species, pdep_network, requires_rms=requires_rms)
 
             else:
                 raise TypeError(
@@ -668,7 +670,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
                         if len(products) == 1 and products[0] == species:
                             new_reactions = network.explore_isomer(species)
 
-                            self.process_new_reactions(new_reactions, species, network)
+                            self.process_new_reactions(new_reactions, species, network, requires_rms=requires_rms)
                             network.update_configurations(self)
                             index = 0
                             break
@@ -693,7 +695,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
                     # Identify a core species which was used to generate the reaction
                     # This is only used to determine the reaction direction for processing
                     spc = spcTuple[0]
-                    self.process_new_reactions(rxnList, spc)
+                    self.process_new_reactions(rxnList, spc, requires_rms=requires_rms)
 
         ################################################################
         # Begin processing the new species and reactions
@@ -705,7 +707,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
 
         # Do thermodynamic filtering
         if not np.isinf(self.thermo_tol_keep_spc_in_edge) and self.new_species_list != []:
-            self.thermo_filter_species(self.new_species_list)
+            self.thermo_filter_species(self.new_species_list, requires_rms=requires_rms)
 
         # Update unimolecular (pressure dependent) reaction networks
         if self.pressure_dependence:
@@ -802,7 +804,7 @@ def clear_surface_adjustments(self):
         self.new_surface_spcs_loss = set()
         self.new_surface_rxns_loss = set()
 
-    def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True):
+    def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True, requires_rms=False):
         """
         Process a list of newly-generated reactions involving the new core
         species or explored isomer `new_species` in network `pdep_network`.
@@ -824,12 +826,12 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g
                     if spec not in self.core.species:
                         all_species_in_core = False
                         if spec not in self.edge.species:
-                            self.add_species_to_edge(spec)
+                            self.add_species_to_edge(spec, requires_rms=requires_rms)
                 for spec in rxn.products:
                     if spec not in self.core.species:
                         all_species_in_core = False
                         if spec not in self.edge.species:
-                            self.add_species_to_edge(spec)
+                            self.add_species_to_edge(spec, requires_rms=requires_rms)
 
             isomer_atoms = sum([len(spec.molecule[0].atoms) for spec in rxn.reactants])
 
@@ -857,9 +859,9 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g
                     # The reaction is not new, so it should already be in the core or edge
                     continue
                 if all_species_in_core:
-                    self.add_reaction_to_core(rxn)
+                    self.add_reaction_to_core(rxn, requires_rms=requires_rms)
                 else:
-                    self.add_reaction_to_edge(rxn)
+                    self.add_reaction_to_edge(rxn, requires_rms=requires_rms)
             else:
                 # Add the reaction to the appropriate unimolecular reaction network
                 # If pdep_network is not None then that will be the network the
@@ -1102,7 +1104,7 @@ def log_enlarge_summary(
         logging.info("    The model edge has {0:d} species and {1:d} reactions".format(edge_species_count, edge_reaction_count))
         logging.info("")
 
-    def add_species_to_core(self, spec):
+    def add_species_to_core(self, spec, requires_rms=False):
         """
         Add a species `spec` to the reaction model core (and remove from edge if
         necessary). This function also moves any reactions in the edge that gain
@@ -1120,7 +1122,7 @@ def add_species_to_core(self, spec):
             if spec in self.edge.species:
                 # remove forbidden species from edge
                 logging.info("Species {0} was Forbidden and not added to Core...Removing from Edge.".format(spec))
-                self.remove_species_from_edge(self.reaction_systems, spec)
+                self.remove_species_from_edge(self.reaction_systems, spec, requires_rms=requires_rms)
                 # remove any empty pdep networks as a result of species removal
                 if self.pressure_dependence:
                     self.remove_empty_pdep_networks()
@@ -1132,7 +1134,8 @@ def add_species_to_core(self, spec):
 
         rxn_list = []
         if spec in self.edge.species:
-            self.edge.phase_system.pass_species(spec.label, self.core.phase_system)
+            if requires_rms:
+                self.edge.phase_system.pass_species(spec.label, self.core.phase_system)
             # If species was in edge, remove it
             logging.debug("Removing species %s from edge.", spec)
             self.edge.species.remove(spec)
@@ -1152,31 +1155,26 @@ def add_species_to_core(self, spec):
 
             # Move any identified reactions to the core
             for rxn in rxn_list:
-                self.add_reaction_to_core(rxn)
+                self.add_reaction_to_core(rxn, requires_rms=requires_rms)
                 logging.debug("Moving reaction from edge to core: %s", rxn)
-        else:
-            if spec.molecule[0].contains_surface_site():
-                self.core.phase_system.phases["Surface"].add_species(spec, edge_phase=self.edge.phase_system.phases["Surface"])
-                self.edge.phase_system.species_dict[spec.label] = spec
-                self.core.phase_system.species_dict[spec.label] = spec
-            else:
-                self.core.phase_system.phases["Default"].add_species(spec, edge_phase=self.edge.phase_system.phases["Default"])
-                self.edge.phase_system.species_dict[spec.label] = spec
-                self.core.phase_system.species_dict[spec.label] = spec
+        elif requires_rms:
+            destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default"
+            self.core.phase_system.phases[destination_phase].add_species(spec, edge_phase=self.edge.phase_system.phases[destination_phase])
+            self.edge.phase_system.species_dict[spec.label] = spec
+            self.core.phase_system.species_dict[spec.label] = spec
 
         return rxn_list
 
-    def add_species_to_edge(self, spec):
+    def add_species_to_edge(self, spec, requires_rms=False):
         """
-        Add a species `spec` to the reaction model edge.
+        Add a species `spec` to the reaction model edge and optionally the RMS phase.
         """
         self.edge.species.append(spec)
-        if spec.molecule[0].contains_surface_site():
-            self.edge.phase_system.phases["Surface"].add_species(spec)
-            self.edge.phase_system.species_dict[spec.label] = spec
-        else:
-            self.edge.phase_system.phases["Default"].add_species(spec)
-            self.edge.phase_system.species_dict[spec.label] = spec
+        if not requires_rms:
+            return
+        destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default"
+        self.edge.phase_system.phases[destination_phase].add_species(spec)
+        self.edge.phase_system.species_dict[spec.label] = spec
 
     def set_thermodynamic_filtering_parameters(
         self, Tmax, thermo_tol_keep_spc_in_edge, min_core_size_for_prune, maximum_edge_species, reaction_systems
@@ -1200,7 +1198,7 @@ def set_thermodynamic_filtering_parameters(
         self.reaction_systems = reaction_systems
         self.maximum_edge_species = maximum_edge_species
 
-    def thermo_filter_species(self, spcs):
+    def thermo_filter_species(self, spcs, requires_rms=False):
         """
         checks Gibbs energy of the species in species against the
         maximum allowed Gibbs energy
@@ -1215,13 +1213,13 @@ def thermo_filter_species(self, spcs):
                     "greater than the thermo_tol_keep_spc_in_edge of "
                     "{3} ".format(spc, G, Gn, self.thermo_tol_keep_spc_in_edge)
                 )
-                self.remove_species_from_edge(self.reaction_systems, spc)
+                self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms)
 
         # Delete any networks that became empty as a result of pruning
         if self.pressure_dependence:
             self.remove_empty_pdep_networks()
 
-    def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0):
+    def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0, requires_rms=False):
         """
         removes species from the edge based on their Gibbs energy until maximum_edge_species
         is reached under the constraint that all removed species are older than
@@ -1263,7 +1261,7 @@ def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_
                 logging.info(
                     "Removing species {0} from edge to meet maximum number of edge species, Gibbs " "number is {1}".format(spc, Gns[rInds[i]])
                 )
-                self.remove_species_from_edge(self.reaction_systems, spc)
+                self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms)
 
             # Delete any networks that became empty as a result of pruning
             if self.pressure_dependence:
@@ -1293,7 +1291,7 @@ def remove_empty_pdep_networks(self):
                     del self.network_dict[source]
                 self.network_list.remove(network)
 
-    def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune):
+    def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune, requires_rms=False):
         """
         Remove species from the model edge based on the simulation results from
         the list of `reaction_systems`.
@@ -1373,7 +1371,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
             for index, spec in species_to_prune[0:prune_due_to_rate_counter]:
                 logging.info("Pruning species %s", spec)
                 logging.debug("    %-56s    %10.4e", spec, max_edge_species_rate_ratios[index])
-                self.remove_species_from_edge(reaction_systems, spec)
+                self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms)
         if len(species_to_prune) - prune_due_to_rate_counter > 0:
             logging.info(
                 "Pruning %d species to obtain an edge size of %d species", len(species_to_prune) - prune_due_to_rate_counter, maximum_edge_species
@@ -1381,7 +1379,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
             for index, spec in species_to_prune[prune_due_to_rate_counter:]:
                 logging.info("Pruning species %s", spec)
                 logging.debug("    %-56s    %10.4e", spec, max_edge_species_rate_ratios[index])
-                self.remove_species_from_edge(reaction_systems, spec)
+                self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms)
 
         # Delete any networks that became empty as a result of pruning
         if self.pressure_dependence:
@@ -1389,7 +1387,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
 
         logging.info("")
 
-    def remove_species_from_edge(self, reaction_systems, spec):
+    def remove_species_from_edge(self, reaction_systems, spec, requires_rms=False):
         """
         Remove species `spec` from the reaction model edge.
         """
@@ -1397,11 +1395,12 @@ def remove_species_from_edge(self, reaction_systems, spec):
         # remove the species
         self.edge.species.remove(spec)
         self.index_species_dict.pop(spec.index)
-        self.edge.phase_system.remove_species(spec)
+        if requires_rms:
+            self.edge.phase_system.remove_species(spec)
 
         # clean up species references in reaction_systems
         for reaction_system in reaction_systems:
-            if not isinstance(reaction_system, Reactor):
+            if not requires_rms or not isinstance(reaction_system, RMSReactor):
                 try:
                     reaction_system.species_index.pop(spec)
                 except KeyError:
@@ -1473,7 +1472,7 @@ def remove_species_from_edge(self, reaction_systems, spec):
             self.species_cache.remove(spec)
             self.species_cache.append(None)
 
-    def add_reaction_to_core(self, rxn):
+    def add_reaction_to_core(self, rxn, requires_rms=False):
         """
         Add a reaction `rxn` to the reaction model core (and remove from edge if
         necessary). This function assumes `rxn` has already been checked to
@@ -1482,7 +1481,8 @@ def add_reaction_to_core(self, rxn):
         """
         if rxn not in self.core.reactions:
             self.core.reactions.append(rxn)
-            if rxn not in self.edge.reactions:
+
+            if requires_rms and rxn not in self.edge.reactions:
                 # If a reaction is not in edge but is going to add to core, it is either a seed mechanism or a newly generated reaction where all reactants and products are already in core
                 # If the reaction is in edge, then the corresponding rms_rxn was moved from edge phase to core phase in pass_species already.
                 rms_species_list = self.core.phase_system.get_rms_species_list()
@@ -1499,7 +1499,7 @@ def add_reaction_to_core(self, rxn):
         if rxn in self.edge.reactions:
             self.edge.reactions.remove(rxn)
 
-    def add_reaction_to_edge(self, rxn):
+    def add_reaction_to_edge(self, rxn, requires_rms=False):
         """
         Add a reaction `rxn` to the reaction model edge. This function assumes
         `rxn` has already been checked to ensure it is supposed to be an edge
@@ -1508,6 +1508,8 @@ def add_reaction_to_edge(self, rxn):
         edge).
         """
         self.edge.reactions.append(rxn)
+        if not requires_rms:
+            return
         rms_species_list = self.edge.phase_system.get_rms_species_list()
         species_names = self.edge.phase_system.get_species_names()
         bits = np.array([spc.molecule[0].contains_surface_site() for spc in rxn.reactants + rxn.products])
@@ -1564,7 +1566,7 @@ def get_stoichiometry_matrix(self):
                     stoichiometry[i, j] = nu
         return stoichiometry.tocsr()
 
-    def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
+    def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=False):
         """
         Add all species and reactions from `seed_mechanism`, a
         :class:`KineticsPrimaryDatabase` object, to the model core. If `react`
@@ -1636,9 +1638,9 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
                 # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                 # We should calculate a pressure-dependent rate for it
                 if len(rxn.reactants) == 1:
-                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
+                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms)
                 else:
-                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])
+                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms)
 
         # Perform species constraints and forbidden species checks
 
@@ -1675,7 +1677,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
                 spec.get_liquid_volumetric_mass_transfer_coefficient_data()
                 spec.get_henry_law_constant_data()
 
-            self.add_species_to_core(spec)
+            self.add_species_to_core(spec, requires_rms=requires_rms)
 
         for rxn in self.new_reaction_list:
             if self.pressure_dependence and rxn.is_unimolecular():
@@ -1686,7 +1688,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
                     submit(spec, self.solvent_name)
 
                 rxn.fix_barrier_height(force_positive=True)
-            self.add_reaction_to_core(rxn)
+            self.add_reaction_to_core(rxn, requires_rms=requires_rms)
 
         # Check we didn't introduce unmarked duplicates
         self.mark_chemkin_duplicates()
@@ -1698,7 +1700,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
             new_edge_reactions=[],
         )
 
-    def add_reaction_library_to_edge(self, reaction_library):
+    def add_reaction_library_to_edge(self, reaction_library, requires_rms=False):
         """
         Add all species and reactions from `reaction_library`, a
         :class:`KineticsPrimaryDatabase` object, to the model edge.
@@ -1763,9 +1765,9 @@ def add_reaction_library_to_edge(self, reaction_library):
                 # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                 # We should calculate a pressure-dependent rate for it
                 if len(rxn.reactants) == 1:
-                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
+                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms)
                 else:
-                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])
+                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms)
 
         # Perform species constraints and forbidden species checks
         for spec in self.new_species_list:
@@ -1802,7 +1804,7 @@ def add_reaction_library_to_edge(self, reaction_library):
                 spec.get_liquid_volumetric_mass_transfer_coefficient_data()
                 spec.get_henry_law_constant_data()
 
-            self.add_species_to_edge(spec)
+            self.add_species_to_edge(spec, requires_rms=requires_rms)
 
         for rxn in self.new_reaction_list:
             if not (
@@ -1817,7 +1819,7 @@ def add_reaction_library_to_edge(self, reaction_library):
                 )
             ):
                 # Don't add to the edge library reactions that were already processed
-                self.add_reaction_to_edge(rxn)
+                self.add_reaction_to_edge(rxn, requires_rms=requires_rms)
 
         if self.save_edge_species:
             from rmgpy.chemkin import mark_duplicate_reaction
diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index 0e1398e51a..89b4104d2b 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -917,7 +917,7 @@ def update(self, reaction_model, pdep_settings):
                                              f'from the {rxn.library} library, and was not added to the model')
                                 break
                         else:
-                            reaction_model.add_reaction_to_core(net_reaction)
+                            reaction_model.add_reaction_to_core(net_reaction, requires_rms=True)
                     else:
                         # Check whether netReaction already exists in the edge as a LibraryReaction
                         for rxn in reaction_model.edge.reactions:
@@ -929,7 +929,7 @@ def update(self, reaction_model, pdep_settings):
                                              f'from the {rxn.library} library, and was not added to the model')
                                 break
                         else:
-                            reaction_model.add_reaction_to_edge(net_reaction)
+                            reaction_model.add_reaction_to_edge(net_reaction, requires_rms=True)
 
                 # Set/update the net reaction kinetics using interpolation model
                 kdata = K[:, :, i, j].copy()
diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactionmechanismsimulator_reactors.py
similarity index 99%
rename from rmgpy/rmg/reactors.py
rename to rmgpy/rmg/reactionmechanismsimulator_reactors.py
index c487601f93..3c748d674b 100644
--- a/rmgpy/rmg/reactors.py
+++ b/rmgpy/rmg/reactionmechanismsimulator_reactors.py
@@ -35,26 +35,6 @@
 import logging
 import itertools
 
-if __debug__:
-    try:
-        from os.path import dirname, abspath, join, exists
-
-        path_rms = dirname(dirname(dirname(abspath(__file__))))
-        from julia.api import Julia
-
-        jl = Julia(sysimage=join(path_rms, "rms.so")) if exists(join(path_rms, "rms.so")) else Julia(compiled_modules=False)
-        from pyrms import rms
-        from diffeqpy import de
-        from julia import Main
-    except Exception as e:
-        import warnings
-
-        warnings.warn("Unable to import Julia dependencies, original error: " + str(e), RuntimeWarning)
-else:
-    from pyrms import rms
-    from diffeqpy import de
-    from julia import Main
-
 from rmgpy.species import Species
 from rmgpy.molecule.fragment import Fragment
 from rmgpy.reaction import Reaction
@@ -71,6 +51,20 @@
 from rmgpy.data.kinetics.family import TemplateReaction
 from rmgpy.data.kinetics.depository import DepositoryReaction
 
+NO_JULIA = False
+try:
+    if __debug__:
+        from os.path import dirname, abspath, join, exists
+        from julia.api import Julia
+        path_rms = dirname(dirname(dirname(abspath(__file__))))
+        jl = Julia(sysimage=join(path_rms, "rms.so")) if exists(join(path_rms, "rms.so")) else Julia(compiled_modules=False)
+    from pyrms import rms
+    from diffeqpy import de
+    from julia import Main
+except Exception as e:
+    logging.info("Unable to import Julia dependencies, original error: " + str(e) + ". RMS features will not be available on this execution.")
+    NO_JULIA = True
+
 
 class PhaseSystem:
     """