diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 75deef0616..a2d3e6d91c 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -150,9 +150,11 @@ jobs:
       # RMS installation and linking to Julia
       - name: Install and link Julia dependencies
         timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours).
+        # JULIA_CONDAPKG_EXE points to the existing conda/mamba to avoid JuliaCall from installing their own. See https://juliapy.github.io/PythonCall.jl/stable/pythoncall/#If-you-already-have-a-Conda-environment.
         run: |
-          python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
-          julia -e 'using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator'
+          juliaup status
+          export JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba
+          julia -e 'ENV["JULIA_CONDAPKG_BACKEND"] = "Current"; using Pkg; Pkg.add(Pkg.PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator'
 
       - name: Install Q2DTor
         run: echo "" | make q2dtor
@@ -169,7 +171,7 @@ jobs:
         run: |
           for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment minimal_surface;
           do
-            if python-jl rmg.py test/regression/"$regr_test"/input.py; then
+            if python rmg.py test/regression/"$regr_test"/input.py; then
               echo "$regr_test" "Executed Successfully"
             else
               echo "$regr_test" "Failed to Execute" | tee -a $GITHUB_STEP_SUMMARY
@@ -258,7 +260,7 @@ jobs:
 
             echo "<details>"
             # Compare the edge and core
-            if python-jl scripts/checkModels.py \
+            if python scripts/checkModels.py \
                 "$regr_test-core" \
                 $REFERENCE/"$regr_test"/chemkin/chem_annotated.inp \
                 $REFERENCE/"$regr_test"/chemkin/species_dictionary.txt \
@@ -275,7 +277,7 @@ jobs:
             cat "$regr_test-core.log"
             echo "</details>"
             echo "<details>"
-            if python-jl scripts/checkModels.py \
+            if python scripts/checkModels.py \
                 "$regr_test-edge" \
                 $REFERENCE/"$regr_test"/chemkin/chem_edge_annotated.inp \
                 $REFERENCE/"$regr_test"/chemkin/species_edge_dictionary.txt \
@@ -296,7 +298,7 @@ jobs:
             if [ -f test/regression/"$regr_test"/regression_input.py ];
             then
               echo "<details>"
-              if python-jl rmgpy/tools/regression.py \
+              if python rmgpy/tools/regression.py \
                 test/regression/"$regr_test"/regression_input.py \
                 $REFERENCE/"$regr_test"/chemkin \
                 test/regression/"$regr_test"/chemkin
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 8ef810f67f..9d8309640c 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -63,8 +63,9 @@ jobs:
       - name: Install and link Julia dependencies
         timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours).
         run: |
-          python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
-          julia -e 'using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator'
+          which julia
+          export JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba
+          julia -e 'ENV["JULIA_CONDAPKG_BACKEND"] = "Current"; using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator'
 
       - name: Checkout gh-pages Branch
         uses: actions/checkout@v2
diff --git a/.gitignore b/.gitignore
index d52d5c680f..946eb94d4a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,7 +47,7 @@ nbproject/*
 .vscode/*
 
 # Unit test files
-.coverage
+.coverage*
 testing/*
 htmlcov/*
 
@@ -117,3 +117,4 @@ examples/**/dictionary.txt
 examples/**/reactions.py
 examples/**/RMG_libraries
 examples/**/species
+examples/**/plots
diff --git a/Dockerfile b/Dockerfile
index a263199c89..d237a659c7 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -13,12 +13,12 @@ RUN ln -snf /bin/bash /bin/sh
 #  - libxrender1 required by RDKit
 RUN apt-get update && \
     apt-get install -y \
-    make \
-    gcc \
-    wget \
-    git \
-    g++ \
-    libxrender1 && \
+        make \
+        gcc \
+        wget \
+        git \
+        g++ \
+        libxrender1 && \
     apt-get autoremove -y && \
     apt-get clean -y
 
@@ -64,16 +64,15 @@ ENV PATH="$RUNNER_CWD/RMG-Py:$PATH"
 # setting this env variable fixes an issue with Julia precompilation on Windows
 ENV JULIA_CPU_TARGET="x86-64,haswell,skylake,broadwell,znver1,znver2,znver3,cascadelake,icelake-client,cooperlake,generic"
 RUN make && \
-    julia -e 'using Pkg; Pkg.add(PackageSpec(name="PyCall",rev="master")); Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator' && \
-    python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
+    julia -e 'ENV["JULIA_CONDAPKG_BACKEND"] = "Current"; using Pkg; Pkg.add(Pkg.PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator'
 
 # RMG-Py should now be installed and ready - trigger precompilation and test run
-RUN python-jl rmg.py examples/rmg/minimal/input.py
+RUN python rmg.py examples/rmg/minimal/input.py
 # delete the results, preserve input.py
 RUN mv examples/rmg/minimal/input.py . && \
     rm -rf examples/rmg/minimal/* && \
     mv input.py examples/rmg/minimal/
 
 # when running this image, open an interactive bash terminal inside the conda environment
-RUN echo "source activate rmg_env" > ~/.bashrc
+RUN echo "source activate rmg_env" >~/.bashrc
 ENTRYPOINT ["/bin/bash", "--login"]
diff --git a/Makefile b/Makefile
index c573fc3b47..5fc01aed95 100644
--- a/Makefile
+++ b/Makefile
@@ -59,16 +59,16 @@ decython:
 	find . -name *.pyc -exec rm -f '{}' \;
 
 test-all:
-	python-jl -m pytest
+	python -m pytest
 
 test test-unittests:
-	python-jl -m pytest -m "not functional and not database"
+	python -m pytest -m "not functional and not database"
 
 test-functional:
-	python-jl -m pytest -m "functional"
+	python -m pytest -m "functional"
 
 test-database:
-	python-jl -m pytest -m "database"
+	python -m pytest -m "database"
 
 eg0: all
 	mkdir -p testing/eg0
diff --git a/documentation/Makefile b/documentation/Makefile
index d42279f698..5c30e66cd7 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/documentation/source/users/rmg/installation/anacondaDeveloper.rst b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
index 4a79244d0c..9eb74479fd 100644
--- a/documentation/source/users/rmg/installation/anacondaDeveloper.rst
+++ b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
@@ -143,7 +143,7 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
 
 #. 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 replace/with/path/to/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.
 
@@ -158,7 +158,7 @@ You might have to edit them slightly to match your exact paths. Specifically,
 you will need ``/opt/miniconda3/envs/rmg_env`` to point to where your conda environment is located.
 
 This configuration will allow you to debug the rms_constant_V example, running through
-python-jl. ::
+python. ::
 
         {
             "name": "Python: rmg.py rms_constant_V",
@@ -166,7 +166,7 @@ python-jl. ::
             "request": "launch",
             "cwd": "${workspaceFolder}/",
             "program": "rmg.py",
-            "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+            "python": "/opt/miniconda3/envs/rmg_env/bin/python",
             "args": [
                 "examples/rmg/rms_constant_V/input.py",
             ],
@@ -185,7 +185,7 @@ Open one of the many test files named ``*Test.py`` in ``test/rmgpy`` before you
             "type": "python",
             "request": "launch",
             "program": "/opt/miniconda3/envs/rmg_env/bin/pytest",
-            "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+            "python": "/opt/miniconda3/envs/rmg_env/bin/python",
             "args": [
                 "--capture=no",
                 "--verbose",
@@ -205,7 +205,7 @@ This configuration will allow you to debug running all the database tests.::
             "type": "python",
             "request": "launch",
             "program": "/opt/miniconda3/envs/rmg_env/bin/pytest",
-            "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+            "python": "/opt/miniconda3/envs/rmg_env/bin/python",
             "args": [
                 "--capture=no",
                 "--verbose",
@@ -226,7 +226,7 @@ This configuration will allow you to use the debugger breakpoints inside unit te
             "request": "launch",
             "program": "${file}",
             "purpose": ["debug-test"],
-            "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+            "python": "/opt/miniconda3/envs/rmg_env/bin/python",
             "console": "integratedTerminal",
             "justMyCode": false,
             "env": {"PYTEST_ADDOPTS": "--no-cov",} // without disabling coverage VS Code doesn't stop at breakpoints while debugging because pytest-cov is using the same technique to access the source code being run
@@ -274,7 +274,7 @@ To configure the rest of the settings, find the ``settings.json`` file in your `
 You can use the following settings to configure the pytest framework::
 
     "python.testing.pytestEnabled": true,
-    "python.testing.pytestPath": "python-jl -m pytest",
+    "python.testing.pytestPath": "python -m pytest",
     "python.testing.pytestArgs": [
         "-p", "julia.pytestplugin",
         "--julia-compiled-modules=no",
diff --git a/documentation/source/users/rmg/running.rst b/documentation/source/users/rmg/running.rst
index 0bb21b80a6..aa28563a89 100755
--- a/documentation/source/users/rmg/running.rst
+++ b/documentation/source/users/rmg/running.rst
@@ -10,7 +10,7 @@ Running a basic RMG job is straightforward, as shown in the example below. Howev
 
 Basic run::
 
-	python-jl rmg.py input.py
+	python rmg.py input.py
 
 .. _inputflags:
 
@@ -19,7 +19,7 @@ Input flags
 
 The options for input flags can be found in ``/RMG-Py/rmgpy/util.py``. Running ::
 
- 	python-jl rmg.py -h
+ 	python rmg.py -h
 
 at the command line will print the documentation from ``util.py``, which is reproduced below for convenience::
 
@@ -63,23 +63,23 @@ Some representative example usages are shown below.
 
 Run by restarting from a seed mechanism::
 
-    python-jl rmg.py -r path/to/seed/ input.py
+    python rmg.py -r path/to/seed/ input.py
 
 Run with CPU time profiling::
 
-    python-jl rmg.py -p input.py
+    python rmg.py -p input.py
 
 Run with multiprocessing for reaction generation and QMTP::
 
-    python-jl rmg.py -n <Max number of processes allowed> input.py 
+    python rmg.py -n <Max number of processes allowed> input.py 
 
 Run with setting a limit on the maximum execution time::
 
-	python-jl rmg.py -t <DD:HH:MM:SS> input.py
+	python rmg.py -t <DD:HH:MM:SS> input.py
 
 Run with setting a limit on the maximum number of iterations::
 
-	python-jl rmg.py -i <Max number of desired iterations> input.py
+	python rmg.py -i <Max number of desired iterations> input.py
 
 
 Details on the multiprocessing implementation
diff --git a/environment.yml b/environment.yml
index 64df68a1c2..8d9871ca36 100644
--- a/environment.yml
+++ b/environment.yml
@@ -16,7 +16,8 @@
 #   made dependency list more explicit (@JacksonBurns).
 # - October 16, 2023 Switched RDKit and descripatastorus to conda-forge,
 #   moved diffeqpy to pip and (temporarily) removed chemprop
-#
+# - May 14, 2024 Removed diffeqpy by switching to call SciMLBase and Sundials using JuliaCall
+
 name: rmg_env
 channels:
   - defaults
@@ -49,8 +50,8 @@ dependencies:
   - conda-forge::rdkit >=2022.09.1
 
 # general-purpose external software tools
-  - conda-forge::julia=1.9.1
-  - conda-forge::pyjulia >=0.6
+  - conda-forge::juliaup
+  - conda-forge::pyjuliacall # for calling julia packages
 
 # Python tools
   - python >=3.7
@@ -88,18 +89,8 @@ dependencies:
   - rmg::pydas >=1.0.3
   - rmg::pydqed >=1.0.3
   - rmg::pyrdl
-  - 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
 
diff --git a/examples/arkane/networks/CH2NH2/input.py b/examples/arkane/networks/CH2NH2_mse/input.py
similarity index 100%
rename from examples/arkane/networks/CH2NH2/input.py
rename to examples/arkane/networks/CH2NH2_mse/input.py
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH3NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH3NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH3NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH3NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH3NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH3NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH3NH_to_CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH_to_CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH3NH_to_CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH_to_CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/H.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/H.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/H.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/H.yml
diff --git a/examples/arkane/networks/acetyl+O2/input.py b/examples/arkane/networks/acetyl+O2_cse/input.py
similarity index 100%
rename from examples/arkane/networks/acetyl+O2/input.py
rename to examples/arkane/networks/acetyl+O2_cse/input.py
diff --git a/examples/arkane/networks/acetyl+O2_rs/input.py b/examples/arkane/networks/acetyl+O2_rs/input.py
new file mode 100644
index 0000000000..84a57961fa
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_rs/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+#   Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+    label = 'acetylperoxy',
+    structure = SMILES('CC(=O)O[O]'),
+    E0 = (-34.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+        HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+        HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+        HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'hydroperoxylvinoxy',
+    structure = SMILES('[CH2]C(=O)OO'),
+    E0 = (-32.4,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+        HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+        HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+        HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'acetyl',
+    structure = SMILES('C[C]=O'),
+    E0 = (0.0,'kcal/mol'),  #(-20.5205,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(43.05,"g/mol")),
+        NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+        HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'oxygen',
+    structure = SMILES('[O][O]'),
+    E0 = (0.0,'kcal/mol'),  #(-5.74557,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(32.00,"g/mol")),
+        LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+    ],
+    spinMultiplicity = 3,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'ketene',
+    structure = SMILES('C=C=O'),
+    E0 = (-6.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(42.0106,"g/mol")),
+        NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'lactone',
+    structure = SMILES('C1OC1(=O)'),
+    E0 = (-30.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(58.0055,"g/mol")),
+        NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+
+
+species(
+    label = 'hydroxyl',
+    structure = SMILES('[OH]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(17.0027,"g/mol")),
+        LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'hydroperoxyl',
+    structure = SMILES('O[O]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(32.9977,"g/mol")),
+        NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'nitrogen',
+    structure = SMILES('N#N'),
+    molecularWeight = (28.04,"g/mol"),
+    collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+    reactive = False
+)
+
+################################################################################
+
+transitionState(
+    label = 'entrance1',
+    E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+    label = 'isom1',
+    E0 = (-5.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit1',
+    E0 = (0.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit2',
+    E0 = (-4.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+        NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+        HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+        HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit3',
+    E0 = (-7.2,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+        HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+    label = 'entrance1',
+    reactants = ['acetyl', 'oxygen'],
+    products = ['acetylperoxy'],
+    transitionState = 'entrance1',
+    kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+    label = 'isom1',
+    reactants = ['acetylperoxy'],
+    products = ['hydroperoxylvinoxy'],
+    transitionState = 'isom1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit1',
+    reactants = ['acetylperoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit2',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit2',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit3',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['lactone', 'hydroxyl'],
+    transitionState = 'exit3',
+    tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+    label = 'acetyl + O2',
+    isomers = [
+        'acetylperoxy',
+        'hydroperoxylvinoxy',
+    ],
+    reactants = [
+        ('acetyl', 'oxygen'),
+        ('ketene', 'hydroperoxyl'),
+	('lactone', 'hydroxyl')
+    ],
+    bathGas = {
+        'nitrogen': 1.0,
+    }
+)
+
+#################################################################################
+
+pressureDependence(
+    'acetyl + O2',
+    Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+    Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+    maximumGrainSize = (1.0,'kcal/mol'),
+    minimumGrainCount = 250,
+    method = 'reservoir state',
+    interpolationModel = ('chebyshev', 6, 4),
+    activeJRotor = True,
+)
diff --git a/examples/arkane/networks/acetyl+O2_sls/input.py b/examples/arkane/networks/acetyl+O2_sls/input.py
new file mode 100644
index 0000000000..4595e288fc
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_sls/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+#   Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+    label = 'acetylperoxy',
+    structure = SMILES('CC(=O)O[O]'),
+    E0 = (-34.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+        HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+        HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+        HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'hydroperoxylvinoxy',
+    structure = SMILES('[CH2]C(=O)OO'),
+    E0 = (-32.4,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+        HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+        HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+        HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'acetyl',
+    structure = SMILES('C[C]=O'),
+    E0 = (0.0,'kcal/mol'),  #(-20.5205,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(43.05,"g/mol")),
+        NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+        HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'oxygen',
+    structure = SMILES('[O][O]'),
+    E0 = (0.0,'kcal/mol'),  #(-5.74557,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(32.00,"g/mol")),
+        LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+    ],
+    spinMultiplicity = 3,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'ketene',
+    structure = SMILES('C=C=O'),
+    E0 = (-6.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(42.0106,"g/mol")),
+        NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'lactone',
+    structure = SMILES('C1OC1(=O)'),
+    E0 = (-30.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(58.0055,"g/mol")),
+        NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+
+
+species(
+    label = 'hydroxyl',
+    structure = SMILES('[OH]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(17.0027,"g/mol")),
+        LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'hydroperoxyl',
+    structure = SMILES('O[O]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(32.9977,"g/mol")),
+        NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'nitrogen',
+    structure = SMILES('N#N'),
+    molecularWeight = (28.04,"g/mol"),
+    collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+    reactive = False
+)
+
+################################################################################
+
+transitionState(
+    label = 'entrance1',
+    E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+    label = 'isom1',
+    E0 = (-5.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit1',
+    E0 = (0.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit2',
+    E0 = (-4.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+        NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+        HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+        HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit3',
+    E0 = (-7.2,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+        HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+    label = 'entrance1',
+    reactants = ['acetyl', 'oxygen'],
+    products = ['acetylperoxy'],
+    transitionState = 'entrance1',
+    kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+    label = 'isom1',
+    reactants = ['acetylperoxy'],
+    products = ['hydroperoxylvinoxy'],
+    transitionState = 'isom1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit1',
+    reactants = ['acetylperoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit2',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit2',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit3',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['lactone', 'hydroxyl'],
+    transitionState = 'exit3',
+    tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+    label = 'acetyl + O2',
+    isomers = [
+        'acetylperoxy',
+        'hydroperoxylvinoxy',
+    ],
+    reactants = [
+        ('acetyl', 'oxygen'),
+        ('ketene', 'hydroperoxyl'),
+	('lactone', 'hydroxyl')
+    ],
+    bathGas = {
+        'nitrogen': 1.0,
+    }
+)
+
+#################################################################################
+
+pressureDependence(
+    'acetyl + O2',
+    Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+    Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+    maximumGrainSize = (1.0,'kcal/mol'),
+    minimumGrainCount = 250,
+    method = 'simulation least squares',
+    interpolationModel = ('chebyshev', 6, 4),
+    activeJRotor = True,
+)
diff --git a/examples/arkane/networks/acetyl+O2_slse/input.py b/examples/arkane/networks/acetyl+O2_slse/input.py
new file mode 100644
index 0000000000..46114439b4
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_slse/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+#   Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+    label = 'acetylperoxy',
+    structure = SMILES('CC(=O)O[O]'),
+    E0 = (-34.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+        HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+        HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+        HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'hydroperoxylvinoxy',
+    structure = SMILES('[CH2]C(=O)OO'),
+    E0 = (-32.4,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+        HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+        HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+        HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    molecularWeight = (75.04,"g/mol"),
+    collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+    energyTransferModel = SingleExponentialDown(
+        alpha0 = (0.5718,'kcal/mol'),
+        T0 = (300,'K'),
+        n = 0.85,
+    ),
+)
+
+species(
+    label = 'acetyl',
+    structure = SMILES('C[C]=O'),
+    E0 = (0.0,'kcal/mol'),  #(-20.5205,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(43.05,"g/mol")),
+        NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+        HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'oxygen',
+    structure = SMILES('[O][O]'),
+    E0 = (0.0,'kcal/mol'),  #(-5.74557,"kJ/mol")
+    modes = [
+        IdealGasTranslation(mass=(32.00,"g/mol")),
+        LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+    ],
+    spinMultiplicity = 3,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'ketene',
+    structure = SMILES('C=C=O'),
+    E0 = (-6.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(42.0106,"g/mol")),
+        NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+        HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'lactone',
+    structure = SMILES('C1OC1(=O)'),
+    E0 = (-30.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(58.0055,"g/mol")),
+        NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+    ],
+    spinMultiplicity = 1,
+    opticalIsomers = 1,
+)
+
+
+
+species(
+    label = 'hydroxyl',
+    structure = SMILES('[OH]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(17.0027,"g/mol")),
+        LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'hydroperoxyl',
+    structure = SMILES('O[O]'),
+    E0 = (0.0,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(32.9977,"g/mol")),
+        NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+)
+
+species(
+    label = 'nitrogen',
+    structure = SMILES('N#N'),
+    molecularWeight = (28.04,"g/mol"),
+    collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+    reactive = False
+)
+
+################################################################################
+
+transitionState(
+    label = 'entrance1',
+    E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+    label = 'isom1',
+    E0 = (-5.8,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit1',
+    E0 = (0.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit2',
+    E0 = (-4.6,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+        NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+        HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+        HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+        HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+    label = 'exit3',
+    E0 = (-7.2,'kcal/mol'),
+    modes = [
+        IdealGasTranslation(mass=(75.04,"g/mol")),
+        NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+        HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+        HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+    ],
+    spinMultiplicity = 2,
+    opticalIsomers = 1,
+    frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+    label = 'entrance1',
+    reactants = ['acetyl', 'oxygen'],
+    products = ['acetylperoxy'],
+    transitionState = 'entrance1',
+    kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+    label = 'isom1',
+    reactants = ['acetylperoxy'],
+    products = ['hydroperoxylvinoxy'],
+    transitionState = 'isom1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit1',
+    reactants = ['acetylperoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit1',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit2',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['ketene', 'hydroperoxyl'],
+    transitionState = 'exit2',
+    tunneling = 'Eckart',
+)
+
+reaction(
+    label = 'exit3',
+    reactants = ['hydroperoxylvinoxy'],
+    products = ['lactone', 'hydroxyl'],
+    transitionState = 'exit3',
+    tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+    label = 'acetyl + O2',
+    isomers = [
+        'acetylperoxy',
+        'hydroperoxylvinoxy',
+    ],
+    reactants = [
+        ('acetyl', 'oxygen'),
+        ('ketene', 'hydroperoxyl'),
+	('lactone', 'hydroxyl')
+    ],
+    bathGas = {
+        'nitrogen': 1.0,
+    }
+)
+
+#################################################################################
+
+pressureDependence(
+    'acetyl + O2',
+    Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+    Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+    maximumGrainSize = (1.0,'kcal/mol'),
+    minimumGrainCount = 250,
+    method = 'simulation least squares eigen',
+    interpolationModel = ('chebyshev', 6, 4),
+    activeJRotor = True,
+)
diff --git a/examples/arkane/networks/n-butanol/input.py b/examples/arkane/networks/n-butanol_msc/input.py
similarity index 100%
rename from examples/arkane/networks/n-butanol/input.py
rename to examples/arkane/networks/n-butanol_msc/input.py
diff --git a/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb b/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
index 700eb3cac3..371bfc2cad 100644
--- a/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
+++ b/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
@@ -49,7 +49,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/catalysis/ch4_o2/input.py"
+    "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/catalysis/ch4_o2/input.py"
    ]
   },
   {
diff --git a/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb b/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
index ba6b3b6a98..44befc780f 100644
--- a/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
+++ b/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
@@ -40,7 +40,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
+    "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
    ]
   },
   {
diff --git a/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb b/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
index 7540b21b01..8ca7eb0029 100644
--- a/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
+++ b/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
@@ -40,7 +40,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
+    "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
    ]
   },
   {
diff --git a/pytest.ini b/pytest.ini
index b17574cc17..7879baabed 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -17,7 +17,6 @@ filterwarnings =
 # -n auto # will use all available cores, if you first pip install pytest-xdist, but it doesn't work well with julia (so add --no-julia? or give up?)
 addopts = 
     --keep-duplicates
-    -p julia.pytestplugin 
     -vv 
     --ignore test/regression
     --cov=arkane --cov=rmgpy
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/rmgpy/pdep/sls.py b/rmgpy/pdep/sls.py
index b183a0bba8..a3474ac6e5 100644
--- a/rmgpy/pdep/sls.py
+++ b/rmgpy/pdep/sls.py
@@ -31,20 +31,17 @@
 Contains functionality for directly simulating the master equation
 and implementing the SLS master equation reduction method
 """
-import sys
 
-from diffeqpy import de
-from julia import Main
-import scipy.sparse as sparse
+from juliacall import Main
+Main.seval("using ReactionMechanismSimulator.SciMLBase")
+Main.seval("using ReactionMechanismSimulator.Sundials")
 import numpy as np
 import scipy.linalg
-import mpmath
 import scipy.optimize as opt
 
 import rmgpy.constants as constants
-from rmgpy.pdep.reaction import calculate_microcanonical_rate_coefficient
 from rmgpy.pdep.me import generate_full_me_matrix, states_to_configurations
-from rmgpy.statmech.translation import IdealGasTranslation
+from rmgpy.rmg.reactors import to_julia
 
 
 def get_initial_condition(network, x0, indices):
@@ -85,30 +82,36 @@ def get_initial_condition(network, x0, indices):
 
 
 def solve_me(M, p0, t):
-    f = Main.eval(
+    f = Main.seval(
         """
-function f(u, M, t)
-    return M*u
+function f(du, u, M, t)
+    du .= M * u
+    return du
 end"""
     )
-    jac = Main.eval(
+    jac = Main.seval(
         """
-function jac(u, M, t)
-    return M
+function jac(J, u, M, t)
+    J .= M
+    return J
 end"""
     )
+    p0 = to_julia(p0)
+    M = to_julia(M)
     tspan = (0.0, t)
-    fcn = de.ODEFunction(f, jac=jac)
-    prob = de.ODEProblem(fcn, p0, tspan, M)
-    sol = de.solve(prob, solver=de.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
+    fcn = Main.ODEFunction(f, jac=jac)
+    prob = Main.ODEProblem(fcn, p0, tspan, M)
+    sol = Main.solve(prob, Main.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
     return sol
 
 
 def solve_me_fcns(f, jac, M, p0, t):
+    p0 = to_julia(p0)
+    M = to_julia(M)
     tspan = (0.0, t)
-    fcn = de.ODEFunction(f, jac=jac)
-    prob = de.ODEProblem(fcn, p0, tspan, M)
-    sol = de.solve(prob, solver=de.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
+    fcn = Main.ODEFunction(f, jac=jac)
+    prob = Main.ODEProblem(fcn, p0, tspan, M)
+    sol = Main.solve(prob, Main.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
     return sol
 
 
@@ -150,16 +153,18 @@ def get_rate_coefficients_SLS(network, T, P, method="mexp", neglect_high_energy_
     tau = np.abs(1.0 / fastest_reaction)
 
     if method == "ode":
-        f = Main.eval(
+        f = Main.seval(
             """
-function f(u,M,t)
-    return M*u
+function f(du, u, M, t)
+    du .= M * u
+    return du
 end"""
         )
-        jac = Main.eval(
+        jac = Main.seval(
             """
-function jac(u,M,t)
-    return M
+function jac(J, u, M, t)
+    J .= M
+    return J
 end"""
         )
 
@@ -247,7 +252,7 @@ def run_equilibrium(network, channel):
     xseq = []
     dxdtseq = []
 
-    # Single domainant source simulations
+    # Single dominant source simulations
     for i, isomer in enumerate(isomers):
         xsout, dxdtout = run_single_source(network, isomer)
         for j in range(len(xsout)):
diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactors.py
index c487601f93..722814a5bc 100644
--- a/rmgpy/rmg/reactors.py
+++ b/rmgpy/rmg/reactors.py
@@ -35,25 +35,12 @@
 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
+import juliacall
+from juliacall import Main
+Main.seval("using PythonCall")
+Main.seval("using Sundials")
+rms = juliacall.newmodule("RMS")
+rms.seval("using ReactionMechanismSimulator")
 
 from rmgpy.species import Species
 from rmgpy.molecule.fragment import Fragment
@@ -72,6 +59,31 @@
 from rmgpy.data.kinetics.depository import DepositoryReaction
 
 
+def to_julia(obj):
+    """
+    Convert native Python object to julia object. If the object is a Python dict, it will be converted to a Julia Dict. If the object is a Python list, it will be converted to a Julia Vector. If the object is a 1-d numpy array, it will be converted to a Julia Array. If the object is a n-d (n > 1) numpy array, it will be converted to a Julia Matrix.
+    Otherwise, the object will be returned as is. Doesn't support RMG objects.
+
+    Parameters
+    ----------
+    obj : dict | list | np.ndarray | object
+        The native Python object to convert
+
+    Returns
+    -------
+    object : Main.Dict | Main.Vector | Main.Matrix | object
+        The Julia object
+    """
+    if isinstance(obj, dict):
+        return Main.PythonCall.pyconvert(Main.Dict, obj)
+    elif isinstance(obj, (list, np.ndarray)):
+        if obj.getattr("shape", False) and len(obj.shape) > 1:
+            return Main.PythonCall.pyconvert(Main.Matrix, obj)
+        return Main.PythonCall.pyconvert(Main.Vector, obj)
+    else: # Other native Python project does not need special conversion.
+        return obj
+
+
 class PhaseSystem:
     """
     Class for tracking and managing all the phases and interfaces of species/reactions
@@ -151,8 +163,11 @@ def pass_species(self, label, phasesys):
 
         rxnlist = []
         for i, rxn in enumerate(self.phases[phase_label].reactions):
-            if (spc.name in [spec.name for spec in rxn.reactants + rxn.products]) and all(
-                [spec.name in phasesys.species_dict for spec in rxn.reactants + rxn.products]
+            reactants = Main.pylist(rxn.reactants)
+            products = Main.pylist(rxn.products)
+            reacs_and_prods = reactants + products
+            if (spc.name in [spec.name for spec in reacs_and_prods]) and all(
+                [spec.name in phasesys.species_dict for spec in reacs_and_prods]
             ):
                 rxnlist.append(rxn)
 
@@ -166,8 +181,8 @@ def pass_species(self, label, phasesys):
             for i, rxn in enumerate(interface.reactions):
                 if (
                     (spc in rxn.reactants or spc in rxn.products)
-                    and all([spec in phasesys.interfaces[key].species for spec in rxn.reactants])
-                    and all([spec in phasesys.interfaces[key].species for spec in rxn.products])
+                    and all(spec.name in phasesys.species_dict for spec in rxn.reactants)
+                    and all(spec.name in phasesys.species_dict for spec in rxn.products)
                 ):
                     rxnlist.append(rxn)
 
@@ -446,7 +461,7 @@ def simulate(self, model_settings, simulator_settings, conditions):
             model_settings.tol_rxn_to_core_deadend_radical,
             atol=simulator_settings.atol,
             rtol=simulator_settings.rtol,
-            solver=de.CVODE_BDF(),
+            solver=Main.Sundials.CVODE_BDF(),
         )
 
         return (
@@ -472,7 +487,7 @@ def generate_reactor(self, phase_system):
         """
         phase = phase_system.phases["Default"]
         ig = rms.IdealGas(phase.species, phase.reactions)
-        domain, y0, p = rms.ConstantVDomain(phase=ig, initialconds=self.initial_conditions)
+        domain, y0, p = rms.ConstantVDomain(phase=ig, initialconds=to_julia(self.initial_conditions))
         react = rms.Reactor(domain, y0, (0.0, self.tf), p=p)
         return react, domain, [], p
 
@@ -492,15 +507,15 @@ def generate_reactor(self, phase_system):
         surf = rms.IdealSurface(surf.species, surf.reactions, surf.site_density, name="surface")
         liq_constant_species = [cspc for cspc in self.const_spc_names if cspc in [spc.name for spc in liq.species]]
         cat_constant_species = [cspc for cspc in self.const_spc_names if cspc in [spc.name for spc in surf.species]]
-        domainliq, y0liq, pliq = rms.ConstantTVDomain(phase=liq, initialconds=self.initial_conditions["liquid"], constantspecies=liq_constant_species)
+        domainliq, y0liq, pliq = rms.ConstantTVDomain(phase=liq, initialconds=to_julia(self.initial_conditions["liquid"]), constantspecies=to_julia(liq_constant_species))
         domaincat, y0cat, pcat = rms.ConstantTAPhiDomain(
-            phase=surf, initialconds=self.initial_conditions["surface"], constantspecies=cat_constant_species
+            phase=surf, initialconds=to_julia(self.initial_conditions["surface"]), constantspecies=to_julia(cat_constant_species),
         )
         if interface.reactions == []:
             inter, pinter = rms.ReactiveInternalInterfaceConstantTPhi(
                 domainliq,
                 domaincat,
-                Main.eval("using ReactionMechanismSimulator; Vector{ElementaryReaction}()"),
+                Main.seval("using ReactionMechanismSimulator; Vector{ElementaryReaction}()"),
                 self.initial_conditions["liquid"]["T"],
                 self.initial_conditions["surface"]["A"],
             )
@@ -536,19 +551,19 @@ def generate_reactor(self, phase_system):
         """
         phase = phase_system.phases["Default"]
         liq = rms.IdealDiluteSolution(phase.species, phase.reactions, phase.solvent)
-        domain, y0, p = rms.ConstantTVDomain(phase=liq, initialconds=self.initial_conditions, constantspecies=self.const_spc_names)
+        domain, y0, p = rms.ConstantTVDomain(phase=liq, initialconds=to_julia(self.initial_conditions), constantspecies=to_julia(self.const_spc_names))
 
         interfaces = []
 
         if self.inlet_conditions:
             inlet_conditions = {key: value for (key, value) in self.inlet_conditions.items() if key != "F"}
             total_molar_flow_rate = self.inlet_conditions["F"]
-            inlet = rms.Inlet(domain, inlet_conditions, Main.eval("x->" + str(total_molar_flow_rate)))
+            inlet = rms.Inlet(domain, to_julia(inlet_conditions), Main.seval("x->" + str(total_molar_flow_rate)))
             interfaces.append(inlet)
 
         if self.outlet_conditions:
             total_volumetric_flow_rate = self.outlet_conditions["Vout"]
-            outlet = rms.VolumetricFlowRateOutlet(domain, Main.eval("x->" + str(total_volumetric_flow_rate)))
+            outlet = rms.VolumetricFlowRateOutlet(domain, Main.seval("x->" + str(total_volumetric_flow_rate)))
             interfaces.append(outlet)
 
         if self.evap_cond_conditions:
@@ -571,7 +586,7 @@ def generate_reactor(self, phase_system):
         """
         phase = phase_system.phases["Default"]
         ig = rms.IdealGas(phase.species, phase.reactions)
-        domain, y0, p = rms.ConstantTPDomain(phase=ig, initialconds=self.initial_conditions)
+        domain, y0, p = rms.ConstantTPDomain(phase=ig, initialconds=to_julia(self.initial_conditions))
         react = rms.Reactor(domain, y0, (0.0, self.tf), p=p)
         return react, domain, [], p
 
@@ -591,21 +606,21 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
         Ea = obj._Ea.value_si
         return rms.Arrhenius(A, n, Ea, rms.EmptyRateUncertainty())
     elif isinstance(obj, PDepArrhenius):
-        Ps = obj._pressures.value_si
-        arrs = [to_rms(arr) for arr in obj.arrhenius]
+        Ps = to_julia(obj._pressures.value_si)
+        arrs = to_julia([to_rms(arr) for arr in obj.arrhenius])
         return rms.PdepArrhenius(Ps, arrs, rms.EmptyRateUncertainty())
     elif isinstance(obj, MultiArrhenius):
-        arrs = [to_rms(arr) for arr in obj.arrhenius]
+        arrs = to_julia([to_rms(arr) for arr in obj.arrhenius])
         return rms.MultiArrhenius(arrs, rms.EmptyRateUncertainty())
     elif isinstance(obj, MultiPDepArrhenius):
-        parrs = [to_rms(parr) for parr in obj.arrhenius]
+        parrs = to_julia([to_rms(parr) for parr in obj.arrhenius])
         return rms.MultiPdepArrhenius(parrs, rms.EmptyRateUncertainty())
     elif isinstance(obj, Chebyshev):
         Tmin = obj.Tmin.value_si
         Tmax = obj.Tmax.value_si
         Pmin = obj.Pmin.value_si
         Pmax = obj.Pmax.value_si
-        coeffs = obj.coeffs.value_si.tolist()
+        coeffs = to_julia(obj.coeffs.value_si)
         return rms.Chebyshev(coeffs, Tmin, Tmax, Pmin, Pmax)
     elif isinstance(obj, ThirdBody):
         arrstr = arrhenius_to_julia_string(obj.arrheniusLow)
@@ -614,7 +629,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
         for key, value in efficiencies.items():
             dstr += '"' + key + '"' "=>" + str(value) + ","
         dstr += "])"
-        return Main.eval(
+        return Main.seval(
             "using ReactionMechanismSimulator; ThirdBody(" + arrstr + ", Dict{Int64,Float64}([]), " + dstr + "," + "EmptyRateUncertainty())"
         )
     elif isinstance(obj, Lindemann):
@@ -625,7 +640,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
         for key, value in efficiencies.items():
             dstr += '"' + key + '"' "=>" + str(value) + ","
         dstr += "])"
-        return Main.eval(
+        return Main.seval(
             "using ReactionMechanismSimulator; Lindemann("
             + arrhigh
             + ","
@@ -648,7 +663,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
         for key, value in efficiencies.items():
             dstr += '"' + key + '"' "=>" + str(value) + ","
         dstr += "])"
-        return Main.eval(
+        return Main.seval(
             "using ReactionMechanismSimulator; Troe("
             + arrhigh
             + ","
@@ -693,6 +708,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
                 atomnums[atm.element.symbol] += 1
             else:
                 atomnums[atm.element.symbol] = 1
+        atomnums = to_julia(atomnums)
         bondnum = len(mol.get_all_edges())
         
         if not obj.molecule[0].contains_surface_site():
@@ -701,12 +717,12 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
             th = obj.get_thermo_data()
             thermo = to_rms(th)
             if obj.henry_law_constant_data:
-                kH = rms.TemperatureDependentHenryLawConstant(Ts=obj.henry_law_constant_data.Ts, kHs=obj.henry_law_constant_data.kHs)
+                kH = rms.TemperatureDependentHenryLawConstant(Ts=to_julia(obj.henry_law_constant_data.Ts), kHs=to_julia(obj.henry_law_constant_data.kHs))
             else:
                 kH = rms.EmptyHenryLawConstant()
             if obj.liquid_volumetric_mass_transfer_coefficient_data:
                 kLA = rms.TemperatureDependentLiquidVolumetricMassTransferCoefficient(
-                    Ts=obj.liquid_volumetric_mass_transfer_coefficient_data.Ts, kLAs=obj.liquid_volumetric_mass_transfer_coefficient_data.kLAs
+                    Ts=to_julia(obj.liquid_volumetric_mass_transfer_coefficient_data.Ts), kLAs=to_julia(obj.liquid_volumetric_mass_transfer_coefficient_data.kLAs)
                 )
             else:
                 kLA = rms.EmptyLiquidVolumetricMassTransferCoefficient()
@@ -748,10 +764,10 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
                 comment=obj.thermo.comment,
             )
     elif isinstance(obj, Reaction):
-        reactantinds = [species_names.index(spc.label) for spc in obj.reactants]
-        productinds = [species_names.index(spc.label) for spc in obj.products]
-        reactants = [rms_species_list[i] for i in reactantinds]
-        products = [rms_species_list[i] for i in productinds]
+        reactantinds = to_julia([species_names.index(spc.label) for spc in obj.reactants])
+        productinds = to_julia([species_names.index(spc.label) for spc in obj.products])
+        reactants = to_julia([rms_species_list[i] for i in reactantinds])
+        products = to_julia([rms_species_list[i] for i in productinds])
         kinetics = to_rms(obj.kinetics, species_names=species_names, rms_species_list=rms_species_list, rmg_species=rmg_species)
         radchange = sum([spc.molecule[0].multiplicity - 1 for spc in obj.products]) - sum([spc.molecule[0].multiplicity - 1 for spc in obj.reactants])
         electronchange = 0  # for now
@@ -765,7 +781,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
             electronchange=electronchange,
             radicalchange=radchange,
             reversible=obj.reversible,
-            pairs=[],
+            pairs=to_julia([]),
             comment=obj.kinetics.comment,
         )
     elif isinstance(obj, SolventData):
diff --git a/test/arkane/arkaneFunctionalInputTest.py b/test/arkane/arkaneFunctionalInputTest.py
index 89f05f6c2e..8e1bcdb1b5 100644
--- a/test/arkane/arkaneFunctionalInputTest.py
+++ b/test/arkane/arkaneFunctionalInputTest.py
@@ -327,7 +327,7 @@ def test_load_input_file(self):
             "examples",
             "arkane",
             "networks",
-            "acetyl+O2",
+            "acetyl+O2_cse",
             "input.py",
         )
         (
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
 
 """