diff --git a/.conda/build.sh b/.conda/build.sh
index c942c713eb..c0ec39927a 100644
--- a/.conda/build.sh
+++ b/.conda/build.sh
@@ -1,6 +1,15 @@
-# Install RMG
+# conda build script
+#
+# conda executes this script during the build process to make the rmg
+# binary
 make install
 
-# lazy "install" of everything in our 'external' folder.
-# most of which should probably be elsewhere
-cp -R ${SRC_DIR}/external ${SP_DIR}
+# The below code does not work because the Julia programming langauge is a very interesting technical demo that is completely unfit for actual
+# deployment in any real world projects.
+# 
+# RMG will be shipped as a pure python package, and then RMS installed by the user.
+#
+# export PYTHON=$PREFIX/bin/python
+# export PYTHONPATH=$SRC_DIR:$PYTHONPATH
+# python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
+# julia -e 'using Pkg; Pkg.add(PackageSpec(name="Functors",version="0.4.3")); Pkg.pin("Functors"); Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="for_rmg")); using ReactionMechanismSimulator'
diff --git a/.conda/conda_build_config.yaml b/.conda/conda_build_config.yaml
index 6a0cca7f0c..4fe9a38c7b 100644
--- a/.conda/conda_build_config.yaml
+++ b/.conda/conda_build_config.yaml
@@ -3,6 +3,3 @@ python:
 numpy:
   - 1.15
 
-# Specifically for Travis build. Should change if building locally.
-CONDA_BUILD_SYSROOT:
-  - /Library/Developer/CommandLineTools/SDKs/MacOSX10.15.sdk # [osx] 
diff --git a/.conda/meta.yaml b/.conda/meta.yaml
index b1df7678e6..cb54c8beda 100644
--- a/.conda/meta.yaml
+++ b/.conda/meta.yaml
@@ -11,66 +11,161 @@ build:
 
 requirements:
   build:
-    - {{ compiler('c') }} # [unix]
+    - {{ compiler('c') }}
   host:
+    - cairo
+    - cairocffi
+    - ffmpeg
+    - xlrd
+    - xlwt
+    - h5py
+    - graphviz
+    - markupsafe
+    - psutil
+    - ncurses
+    - suitesparse
+
+  # external software tools for chemistry
+    - coolprop
+    - cantera=2.6
+    - mopac
+    - cclib >=1.6.3,!=1.8.0
+    - openbabel >=3
+
+  # Python tools
+    - python >=3.7
+    - coverage
     - cython >=0.25.2
+    - scikit-learn
+    - scipy
+    - numpy >=1.10.0
+    - pydot
+    - jinja2
+    - jupyter
+    - pymongo
+    - pyparsing
+    - pyyaml
+    - networkx
+    - matplotlib >=1.5
+    - mpmath
+    - pandas
+
+  # packages we maintain
+    - rmgdatabase
+    - gprof2dot
     - lpsolve55
-    - numpy
-    - openbabel >=3
-    - pydas >=1.0.2
-    - pydqed >=1.0.1
+    - muq2
+    - numdifftools
+    - pydas >=1.0.3
+    - pydqed >=1.0.3
     - pyrdl
-    - python
     - quantities
-    - rdkit >=2018
-    - scipy
-    - setuptools
+    - symmetry
+    - chemprop==0.0.1
+    - rdkit >=2020.03.3.0
   run:
     - cairo
     - cairocffi
-    - cantera >=2.3.0
-    - cclib >=1.6.3
-    - chemprop
+    - ffmpeg
+    - xlrd
+    - xlwt
+    - h5py
+    - graphviz
+    - markupsafe
+    - psutil
+    - ncurses
+    - suitesparse
+
+  # external software tools for chemistry
     - coolprop
+    - cantera=2.6
+    - mopac
+    - cclib >=1.6.3,!=1.8.0
+    - openbabel >=3
+
+  # Python tools
+    - python >=3.7
     - coverage
     - cython >=0.25.2
-    - ffmpeg
-    - gprof2dot
-    - graphviz
-    - h5py
+    - scikit-learn
+    - scipy
+    - numpy >=1.10.0
+    - pydot
     - jinja2
     - jupyter
-    - lpsolve55
-    - markupsafe
+    - pymongo
+    - pyparsing
+    - pyyaml
+    - networkx
     - matplotlib >=1.5
-    - mopac
     - mpmath
+    - pandas
+
+  # packages we maintain
+    - rmgdatabase
+    - gprof2dot
+    - lpsolve55
     - muq2
-    - networkx
-    - nose
     - numdifftools
-    - {{ pin_compatible('numpy') }}
-    - openbabel >=3
-    - pandas
+    - pydas >=1.0.3
+    - pydqed >=1.0.3
+    - pyrdl
+    - quantities
+    - symmetry
+    - chemprop==0.0.1
+    - rdkit >=2020.03.3.0
+test:
+  requires:
+    - cairo
+    - cairocffi
+    - ffmpeg
+    - xlrd
+    - xlwt
+    - h5py
+    - graphviz
+    - markupsafe
     - psutil
-    - pydas >=1.0.2
+    - ncurses
+    - suitesparse
+
+  # external software tools for chemistry
+    - coolprop
+    - cantera=2.6
+    - mopac
+    - cclib >=1.6.3,!=1.8.0
+    - openbabel >=3
+
+  # Python tools
+    - python >=3.7
+    - coverage
+    - cython >=0.25.2
+    - scikit-learn
+    - scipy
+    - numpy >=1.10.0
     - pydot
-    - pydqed >=1.0.1
+    - jinja2
+    - jupyter
     - pymongo
     - pyparsing
-    - pyrdl
-    - python
     - pyyaml
-    - pyzmq
+    - networkx
+    - matplotlib >=1.5
+    - mpmath
+    - pandas
+
+  # packages we maintain
+    - rmgdatabase
+    - gprof2dot
+    - lpsolve55
+    - muq2
+    - numdifftools
+    - pydas >=1.0.3
+    - pydqed >=1.0.3
+    - pyrdl
     - quantities
-    - rdkit >=2018
-    - rmgdatabase >=3.2.0
-    - scikit-learn
-    - scipy
     - symmetry
-    - xlrd
-    - xlwt
-test:
+    - chemprop==0.0.1
+    - rdkit >=2020.03.3.0
   source_files:
     - 'examples/rmg/superminimal'
     - 'examples/arkane/networks/n-butanol'
@@ -78,10 +173,8 @@ test:
     - rmgpy
     - arkane
   commands:
-    - rmg.py examples/rmg/superminimal/input.py # [unix]
-    - Arkane.py examples/arkane/networks/n-butanol/input.py # [unix]
-    - python %SCRIPTS%\rmg.py examples\rmg\superminimal\input.py # [win]
-    - python %SCRIPTS\Arkane.py examples\arkane\networks\n-butanol\input.py # [win]
+    - rmg.py examples/rmg/superminimal/input.py
+    - Arkane.py examples/arkane/networks/n-butanol/input.py
 
 about:
   home: https://github.com/ReactionMechanismGenerator/RMG-Py
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index fa18af19c9..12539a716b 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: |
+          conda install -c conda-forge --override-channels julia>=1.8.5,!=1.9.0 pyjulia>=0.6 
+          conda install -c rmg --override-channels numdifftools pyrms diffeqpy
+
       # list the environment for debugging purposes
       - name: mamba info
         run: |
diff --git a/.github/workflows/conda_build.yml b/.github/workflows/conda_build.yml
deleted file mode 100644
index 8ddbef0877..0000000000
--- a/.github/workflows/conda_build.yml
+++ /dev/null
@@ -1,58 +0,0 @@
-name: Conda Build
-
-on:
-  push:
-    branches:
-      - stable
-jobs:
-  build-linux:
-    runs-on: ubuntu-latest
-    defaults:
-      run:
-        shell: bash -l {0}
-    steps:
-      - uses: actions/checkout@v2
-      - uses: conda-incubator/setup-miniconda@v2
-        with:
-          environment-file: environment.yml
-          python-version: 3.7
-          activate-environment: rmg_env
-      - name: Conda info
-        run: |
-          conda info
-          conda list
-      - name: Build Binary
-        env: 
-          CONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }}
-        run: |
-          conda install -y conda-build 
-          conda install -y anaconda-client 
-          conda config --add channels rmg
-          conda config --set anaconda_upload yes
-          conda build --token $CONDA_TOKEN --user rmg .conda
-  build-osx:
-    runs-on: macos-latest
-    defaults:
-      run:
-        shell: bash -l {0}
-    steps:
-      - uses: actions/checkout@v2
-      - uses: conda-incubator/setup-miniconda@v2
-        with:
-          environment-file: environment.yml
-          python-version: 3.7
-          activate-environment: rmg_env
-      - name: Conda info
-        run: |
-          conda info
-          conda list
-      - name: Build Binary
-        env: 
-          CONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }}
-        run: |
-            conda install -y conda-build 
-            conda install -y anaconda-client 
-            conda config --add channels rmg
-            conda config --set anaconda_upload yes
-            xcrun --show-sdk-path
-            conda build --token $CONDA_TOKEN --user rmg .conda
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 82c694d42b..2a38d443c9 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -31,11 +31,6 @@ jobs:
         run: |
           mamba info
           mamba list
-      - name: Install and link Julia dependencies 
-        run: |
-          julia -e "using Pkg; Pkg.add(PackageSpec(url=\"https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl\", rev=\"main\"))"
-          julia -e "using Pkg; Pkg.add(\"PyCall\"); Pkg.add(\"DifferentialEquations\")"
-          python -c "import julia; julia.install()"
       - name: Install and compile RMG
         run: |
           cd ..
diff --git a/.rmgmolecule_conda/build.sh b/.rmgmolecule_conda/build.sh
new file mode 100644
index 0000000000..b68ece3906
--- /dev/null
+++ b/.rmgmolecule_conda/build.sh
@@ -0,0 +1,2 @@
+# Install rmgmolecule
+python rmgmolecule_setup.py install
diff --git a/.rmgmolecule_conda/conda_build_config.yaml b/.rmgmolecule_conda/conda_build_config.yaml
new file mode 100644
index 0000000000..83e1bc711e
--- /dev/null
+++ b/.rmgmolecule_conda/conda_build_config.yaml
@@ -0,0 +1,2 @@
+numpy:
+  - 1.19
\ No newline at end of file
diff --git a/.rmgmolecule_conda/meta.yaml b/.rmgmolecule_conda/meta.yaml
new file mode 100644
index 0000000000..ce1e04118b
--- /dev/null
+++ b/.rmgmolecule_conda/meta.yaml
@@ -0,0 +1,48 @@
+# meta.yaml - package specification for rmgmolecule subpackage
+package:
+  name: rmgmolecule
+  version: 1.0.0
+
+source:
+  path: ../
+
+build:
+  number: {{ environ.get('GIT_DESCRIBE_NUMBER', 0) }}
+
+requirements:
+  build:
+    - {{ compiler('c') }}
+  host:
+    - cairo
+    - cairocffi
+    - cython >=0.25.2
+    - lpsolve55
+    - numpy
+    - openbabel >=3
+    - pydot
+    - pyrdl
+    - python==3.7
+    - quantities
+    - rdkit >=2018
+  run:
+    - cairo
+    - cairocffi
+    - cython >=0.25.2
+    - lpsolve55
+    - numpy
+    - openbabel >=3
+    - pydot
+    - pyrdl
+    - python==3.7
+    - quantities
+    - rdkit >=2018
+test:
+  imports:
+    - rmgpy.molecule
+  commands:
+    - python -c 'from rmgpy.molecule import Molecule; mol=Molecule().from_smiles("CC")'
+
+about:
+  home: https://github.com/ReactionMechanismGenerator/RMG-Py
+  license: MIT
+  summary: "The ReactionMechanismGenerator Molecule Subpackage"
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
index 803374c02a..149e8b89dd 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -45,6 +45,8 @@ RUN git clone --single-branch --branch main --depth 1 https://github.com/Reactio
 # build the conda environment
 WORKDIR /rmg/RMG-Py
 RUN conda env create --file environment.yml && \
+    conda install -c conda-forge julia>=1.8.5,!=1.9.0 pyjulia>=0.6 && \
+    conda install -c rmg numdifftools pyrms diffeqpy && \
     conda clean --all --yes
 
 # This runs all subsequent commands inside the rmg_env conda environment
diff --git a/documentation/Makefile b/documentation/Makefile
index fb901249e9..cac86744cf 100644
--- a/documentation/Makefile
+++ b/documentation/Makefile
@@ -3,7 +3,7 @@
 
 # You can set these variables from the command line.
 SPHINXOPTS    =
-SPHINXBUILD   = python-jl $$(which sphinx-build)
+SPHINXBUILD   = python $$(which sphinx-build)
 PAPER         =
 BUILDDIR      = build
 
diff --git a/environment.yml b/environment.yml
index 5155a74b09..be687f85b2 100644
--- a/environment.yml
+++ b/environment.yml
@@ -45,10 +45,6 @@ dependencies:
   - conda-forge::cclib >=1.6.3,!=1.8.0
   - conda-forge::openbabel >= 3
 
-# general-purpose external software tools
-  - conda-forge::julia>=1.8.5,!=1.9.0
-  - conda-forge::pyjulia >=0.6
-
 # Python tools
   - python >=3.7
   - coverage
@@ -75,21 +71,20 @@ dependencies:
   - matplotlib >=1.5
   - mpmath
   - pandas
+  - conda-forge::numdifftools
 
 # packages we maintain
   - rmg::gprof2dot
   - rmg::lpsolve55
   - rmg::muq2
-  - rmg::numdifftools
   - rmg::pydas >=1.0.3
   - rmg::pydqed >=1.0.3
   - rmg::pyrdl
-  - rmg::pyrms
   - rmg::quantities
   - rmg::symmetry
 
 # packages we would like to stop maintaining (and why)
-  - rmg::diffeqpy
+  #- 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:
diff --git a/rmg.py b/rmg.py
index 26c27c00a4..0b27098b81 100644
--- a/rmg.py
+++ b/rmg.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python-jl
+#!/usr/bin/env python
 
 ###############################################################################
 #                                                                             #
diff --git a/rmgmolecule_setup.py b/rmgmolecule_setup.py
new file mode 100644
index 0000000000..8614e22ede
--- /dev/null
+++ b/rmgmolecule_setup.py
@@ -0,0 +1,176 @@
+import sys
+import os
+from collections import OrderedDict
+
+try:
+    from distutils.core import setup
+    from distutils.extension import Extension
+except ImportError:
+    print("The distutils package is required to build or install RMG Py.")
+    raise
+
+try:
+    from Cython.Build import cythonize
+    from Cython.Compiler import Options
+except ImportError:
+    print("Cython (http://www.cython.org/) is required to build or install RMG Py.")
+    raise
+
+try:
+    import numpy
+except ImportError:
+    print("NumPy (http://numpy.scipy.org/) is required to build or install RMG Py.")
+    raise
+
+# Create annotated HTML files for each of the Cython modules
+Options.annotate = True
+
+directives = {
+    # Set input language version to python 3
+    "language_level": 3,
+    # Turn on profiling capacity for all Cython modules
+    # 'profile': True,
+    # Embed call signatures in cythonized files - enable when building documentation
+    # 'embedsignature': True,
+}
+
+
+main_ext_modules = [
+    # Molecules and molecular representations
+    Extension(
+        "rmgpy.molecule.atomtype",
+        ["rmgpy/molecule/atomtype.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.element",
+        ["rmgpy/molecule/element.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.graph",
+        ["rmgpy/molecule/graph.pyx"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.group",
+        ["rmgpy/molecule/group.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.molecule",
+        ["rmgpy/molecule/molecule.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.symmetry",
+        ["rmgpy/molecule/symmetry.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.vf2",
+        ["rmgpy/molecule/vf2.pyx"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.converter",
+        ["rmgpy/molecule/converter.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.translator",
+        ["rmgpy/molecule/translator.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.util",
+        ["rmgpy/molecule/util.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.inchi",
+        ["rmgpy/molecule/inchi.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.resonance",
+        ["rmgpy/molecule/resonance.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.pathfinder",
+        ["rmgpy/molecule/pathfinder.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.molecule.kekulize",
+        ["rmgpy/molecule/kekulize.pyx"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.constants",
+        ["rmgpy/constants.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.quantity",
+        ["rmgpy/quantity.py"],
+        include_dirs=["."],
+    ),
+    Extension(
+        "rmgpy.rmgobject",
+        ["rmgpy/rmgobject.pyx"],
+    ),
+]
+
+ext_modules = []
+if "install" in sys.argv:
+    # This is so users can still do simply `python setup.py install`
+    ext_modules.extend(main_ext_modules)
+if "main" in sys.argv:
+    # This is for `python setup.py build_ext main`
+    sys.argv.remove("main")
+    ext_modules.extend(main_ext_modules)
+if "minimal" in sys.argv:
+    # This starts with the full install list, but removes anything that has a pure python mode
+    # i.e. in only includes things whose source is .pyx
+    sys.argv.remove("minimal")
+    temporary_list = []
+    temporary_list.extend(main_ext_modules)
+    for module in temporary_list:
+        for source in module.sources:
+            if os.path.splitext(source)[1] == ".pyx":
+                ext_modules.append(module)
+
+# Remove duplicates while preserving order:
+ext_modules = list(OrderedDict.fromkeys(ext_modules))
+
+scripts = []
+
+modules = ["rmgpy.exceptions", "rmgpy.version"]
+
+__version__ = "1.0.0"
+
+import logging
+
+logging.error(ext_modules)
+# Initiate the build and/or installation
+setup(
+    name="RMG-Py",
+    version=__version__,
+    description="Reaction Mechanism Generator",
+    author="William H. Green and the RMG Team",
+    author_email="rmg_dev@mit.edu",
+    url="http://reactionmechanismgenerator.github.io",
+    packages=[
+        "rmgpy.molecule",
+    ],
+    py_modules=modules,
+    scripts=scripts,
+    ext_modules=cythonize(
+        ext_modules,
+        build_dir="build",
+        compiler_directives=directives,
+    ),
+    include_dirs=[".", numpy.get_include()],
+)
diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py
index 6dd6c8b46b..c79fc19ca8 100644
--- a/rmgpy/molecule/draw.py
+++ b/rmgpy/molecule/draw.py
@@ -61,7 +61,12 @@
 from rdkit.Chem import AllChem
 
 from rmgpy.molecule.molecule import Atom, Molecule
-from rmgpy.qm.molecule import Geometry
+
+Geometry = None
+try:
+    from rmgpy.qm.molecule import Geometry
+except ImportError:
+    logging.info("Unable to import Geometry rmgpy.qm.molecule - feature disabled.")
 
 
 ################################################################################
@@ -437,6 +442,11 @@ def _generate_coordinates(self):
             # Generate the RDkit molecule from the RDkit molecule, use geometry
             # in order to match the atoms in the rdmol with the atoms in the
             # RMG molecule (which is required to extract coordinates).
+            if Geometry is None:
+                raise RuntimeError("""
+Missing rmgpy.qm.molecule.Geometry required for 2D coordinate generation.
+Please install the full version of RMG to use this function.
+                    """)
             self.geometry = Geometry(None, None, self.molecule, None)
 
             rdmol, rd_atom_idx = self.geometry.rd_build()
diff --git a/rmgpy/molecule/fragment.py b/rmgpy/molecule/fragment.py
index 8b30c1e32e..06e25f7346 100644
--- a/rmgpy/molecule/fragment.py
+++ b/rmgpy/molecule/fragment.py
@@ -1020,7 +1020,13 @@ def assign_representative_molecule(self):
 
     def assign_representative_species(self):
 
-        from rmgpy.species import Species
+        try:
+            from rmgpy.species import Species
+        except ImportError as ie:
+            raise RuntimeError(
+                "Cannot call Fragment.assign_representative_species() in rmgmolecule installation."
+                " Please install the full version of RMG."
+                ) from ie
         self.assign_representative_molecule()
         self.species_repr = Species(molecule=[self.mol_repr])
         self.symmetry_number = self.get_symmetry_number()
diff --git a/rmgpy/molecule/inchi.py b/rmgpy/molecule/inchi.py
index 59e49d8522..789a00e0b0 100644
--- a/rmgpy/molecule/inchi.py
+++ b/rmgpy/molecule/inchi.py
@@ -702,7 +702,13 @@ def _convert_3_atom_2_bond_path(start, mol):
     with a number of actions that reflect the changes in bond orders and unpaired
     electrons that the molecule should undergo.
     """
-    from rmgpy.data.kinetics.family import ReactionRecipe
+    try:
+        from rmgpy.data.kinetics.family import ReactionRecipe
+    except ImportError as ie:
+        raise RuntimeError(
+            "Cannot call inchi._convert_3_atom_2_bond_path() in rmgmolecule installation."
+            " Please install the full version of RMG."
+            ) from ie 
 
     def is_valid(mol):
         """Check if total bond order of oxygen atoms is smaller than 4."""
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index e04cd82e2e..f776d63027 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,8 @@ 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.")
         logging.exception(e)
         raise
     finally:
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index 8291db4a13..6fe384fbb6 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -79,7 +79,8 @@
 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
+from rmgpy.rmg.reactionmechanismsimulator_reactors import NO_JULIA
 
 ################################################################################
 
@@ -506,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
@@ -549,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 not NO_JULIA and isinstance(reaction_system, RMSReactor):
                 reaction_system.finish_termination_criteria()
 
         # Load restart seed mechanism (if specified)
@@ -617,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")]:
@@ -694,35 +701,35 @@ 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
+                if ((NO_JULIA or 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 not NO_JULIA and requires_rms:
+            self.attach(RMSWriter(self.output_directory))
 
         if self.generate_output_html:
             self.attach(OutputHTMLWriter(self.output_directory))
@@ -734,7 +741,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 not NO_JULIA and 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))
@@ -748,10 +755,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
 
@@ -778,7 +785,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 not NO_JULIA and 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),
@@ -821,6 +828,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):
@@ -833,7 +841,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")
 
@@ -899,7 +907,7 @@ def execute(self, initialize=True, **kwargs):
                             prune = False
 
                         try:
-                            if isinstance(reaction_system, Reactor):
+                            if not NO_JULIA and requires_rms and isinstance(reaction_system, RMSReactor):
                                 (
                                     terminated,
                                     resurrected,
@@ -992,7 +1000,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
@@ -1001,7 +1009,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 not NO_JULIA and requires_rms and isinstance(reaction_system, RMSReactor):
                                         (
                                             terminated,
                                             resurrected,
@@ -1070,7 +1078,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("")
@@ -1093,13 +1101,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
 
@@ -1126,6 +1135,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()
@@ -1890,7 +1900,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
         """
@@ -1923,6 +1933,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(
@@ -2217,7 +2228,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 not NO_JULIA and 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 951b670d00..e7a1dc2bad 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
@@ -635,13 +637,13 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
                     display(new_species)  # if running in IPython --pylab mode, draws the picture!
 
                 # Add new 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)
 
             elif isinstance(new_object, tuple) and isinstance(new_object[0], PDepNetwork) and self.pressure_dependence:
                 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(
@@ -665,7 +667,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
@@ -690,7 +692,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
@@ -702,7 +704,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:
@@ -799,7 +801,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`.
@@ -821,12 +823,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])
 
@@ -854,9 +856,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
@@ -1111,7 +1113,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
@@ -1129,7 +1131,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()
@@ -1141,7 +1143,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 not NO_JULIA and 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)
@@ -1161,31 +1164,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 not NO_JULIA and 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 NO_JULIA or 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
@@ -1209,7 +1207,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
@@ -1224,13 +1222,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
@@ -1272,7 +1270,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:
@@ -1302,7 +1300,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`.
@@ -1382,7 +1380,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
@@ -1390,7 +1388,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:
@@ -1398,7 +1396,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.
         """
@@ -1406,11 +1404,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 not NO_JULIA and 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 NO_JULIA or not requires_rms or not isinstance(reaction_system, RMSReactor):
                 try:
                     reaction_system.species_index.pop(spec)
                 except KeyError:
@@ -1482,7 +1481,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
@@ -1491,7 +1490,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 not NO_JULIA and 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()
@@ -1508,7 +1508,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
@@ -1517,6 +1517,8 @@ def add_reaction_to_edge(self, rxn):
         edge).
         """
         self.edge.reactions.append(rxn)
+        if NO_JULIA or 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])
@@ -1573,7 +1575,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`
@@ -1641,9 +1643,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
 
@@ -1680,7 +1682,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():
@@ -1691,7 +1693,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()
@@ -1703,7 +1705,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 cc30b86116..d778aabca2 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.reaction import Reaction
 from rmgpy.thermo.nasa import NASAPolynomial, NASA
@@ -70,6 +50,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:
     """
diff --git a/setup.py b/setup.py
index 604a856eb1..17fdab4a88 100644
--- a/setup.py
+++ b/setup.py
@@ -215,7 +215,6 @@
     'scripts/standardizeModelSpeciesNames.py',
     'scripts/thermoEstimator.py',
     'scripts/isotopes.py',
-    'testing/databaseTest.py',
 ]
 
 modules = []
diff --git a/utilities.py b/utilities.py
index c8e7425818..f84963a632 100644
--- a/utilities.py
+++ b/utilities.py
@@ -310,7 +310,7 @@ def update_headers():
     start of each file, be sure to double-check the results to make sure
     important lines aren't accidentally overwritten.
     """
-    shebang = """#!/usr/bin/env python-jl
+    shebang = """#!/usr/bin/env python
 
 """