diff --git a/build-python/fix_rpaths.py b/build-python/fix_rpaths.py
new file mode 100644
index 00000000..34b4f3a8
--- /dev/null
+++ b/build-python/fix_rpaths.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+import glob
+import subprocess
+from collections import OrderedDict
+
+def get_rpaths(binary_path):
+    """Uses otool to extract a list of all LC_RPATH entries."""
+    print(f"--- Checking rpaths for: {binary_path}")
+    rpaths = []
+    try:
+        proc = subprocess.run(
+            ['otool', '-l', binary_path],
+            capture_output=True,
+            text=True,
+            check=True
+        )
+
+        lines = proc.stdout.splitlines()
+        for i, line in enumerate(lines):
+            if "cmd LC_RPATH" in line.strip():
+                if i + 2 < len(lines):
+                    path_line = lines[i + 2].strip()
+                    if path_line.startswith("path "):
+                        # Extract the path, e.g., "path /foo/bar (offset 12)"
+                        rpath = path_line.split(" ")[1]
+                        rpaths.append(rpath)
+
+    except (subprocess.CalledProcessError, FileNotFoundError) as e:
+        print(f"Error running otool: {e}")
+        return []
+
+    return rpaths
+
+def fix_rpaths(binary_path):
+    all_rpaths = get_rpaths(binary_path)
+    if not all_rpaths:
+        print("--- No rpaths found or otool failed.")
+        return
+
+    unique_rpaths = list(OrderedDict.fromkeys(all_rpaths))
+    for rpath in unique_rpaths:
+        print(f" - RPATH: {rpath}")
+
+    if len(all_rpaths) == len(unique_rpaths):
+        print("--- No duplicate rpaths found. Nothing to do.")
+        return
+
+    print(f"--- Found {len(all_rpaths)} rpaths; {len(unique_rpaths)} are unique.")
+    print(f"--- Fixing duplicates in: {binary_path}")
+
+    try:
+        for rpath in all_rpaths:
+            subprocess.run(
+                ['install_name_tool', '-delete_rpath', rpath, binary_path],
+                check=True,
+                capture_output=True
+            )
+
+        for rpath in unique_rpaths:
+            subprocess.run(
+                ['install_name_tool', '-add_rpath', rpath, binary_path],
+                check=True,
+                capture_output=True
+            )
+
+        print("--- Successfully fixed rpaths.")
+
+    except (subprocess.CalledProcessError, FileNotFoundError) as e:
+        print(f"--- Error running install_name_tool: {e}")
+        if e.stderr:
+            print(f"STDERR: {e.stderr.decode()}")
+        if e.stdout:
+            print(f"STDOUT: {e.stdout.decode()}")
+        sys.exit(1) # Fail the install if we can't fix it
+
+def main():
+    if len(sys.argv) != 2:
+        print(f"--- Error: Expected one argument (path to .so file), got {sys.argv}", file=sys.stderr)
+        sys.exit(1)
+
+    # Get the file path directly from the command line argument
+    so_file_path = sys.argv[1]
+
+    if not os.path.exists(so_file_path):
+        print(f"--- Error: File not found at {so_file_path}", file=sys.stderr)
+        sys.exit(1)
+
+    print(f"--- Fixing rpaths for built file: {so_file_path}")
+    fix_rpaths(so_file_path)
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/build-python/meson.build b/build-python/meson.build
index 1da2669f..671ed52c 100644
--- a/build-python/meson.build
+++ b/build-python/meson.build
@@ -14,20 +14,20 @@ py_mod = py_installation.extension_module(
     sources: [
         meson.project_source_root() + '/src/python/bindings.cpp',
         meson.project_source_root() + '/src/python/types/bindings.cpp',
-#        meson.project_source_root() + '/src/python/expectations/bindings.cpp',
-#        meson.project_source_root() + '/src/python/partition/bindings.cpp',
-#        meson.project_source_root() + '/src/python/partition/trampoline/py_partition.cpp',
-#        meson.project_source_root() + '/src/python/reaction/bindings.cpp',
-#        meson.project_source_root() + '/src/python/screening/bindings.cpp',
-#        meson.project_source_root() + '/src/python/screening/trampoline/py_screening.cpp',
-#        meson.project_source_root() + '/src/python/io/bindings.cpp',
-#        meson.project_source_root() + '/src/python/io/trampoline/py_io.cpp',
-#        meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
-#        meson.project_source_root() + '/src/python/engine/bindings.cpp',
-#        meson.project_source_root() + '/src/python/engine/trampoline/py_engine.cpp',
-#        meson.project_source_root() + '/src/python/solver/bindings.cpp',
-#        meson.project_source_root() + '/src/python/solver/trampoline/py_solver.cpp',
-#        meson.project_source_root() + '/src/python/utils/bindings.cpp',
+        meson.project_source_root() + '/src/python/expectations/bindings.cpp',
+        meson.project_source_root() + '/src/python/partition/bindings.cpp',
+        meson.project_source_root() + '/src/python/partition/trampoline/py_partition.cpp',
+        meson.project_source_root() + '/src/python/reaction/bindings.cpp',
+        meson.project_source_root() + '/src/python/screening/bindings.cpp',
+        meson.project_source_root() + '/src/python/screening/trampoline/py_screening.cpp',
+        meson.project_source_root() + '/src/python/io/bindings.cpp',
+        meson.project_source_root() + '/src/python/io/trampoline/py_io.cpp',
+        meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
+        meson.project_source_root() + '/src/python/engine/bindings.cpp',
+        meson.project_source_root() + '/src/python/engine/trampoline/py_engine.cpp',
+        meson.project_source_root() + '/src/python/solver/bindings.cpp',
+        meson.project_source_root() + '/src/python/solver/trampoline/py_solver.cpp',
+        meson.project_source_root() + '/src/python/utils/bindings.cpp',
     ],
     dependencies : gridfire_py_deps,
     cpp_args : ['-UNDEBUG'], # Example: Ensure assertions are enabled if needed
diff --git a/pip_install_mac_patch.sh b/pip_install_mac_patch.sh
new file mode 100755
index 00000000..4f5b7508
--- /dev/null
+++ b/pip_install_mac_patch.sh
@@ -0,0 +1,143 @@
+#!/bin/bash
+# pip_install_mac.sh - Temporary workaround for meson-python duplicate RPATH bug on macOS
+
+set -e
+
+# Color codes for output
+RED='\033[0;31m'
+YELLOW='\033[1;33m'
+GREEN='\033[0;32m'
+NC='\033[0m' # No Color
+
+echo -e "${YELLOW}"
+echo "========================================================================="
+echo "                    TEMPORARY INSTALLATION WORKAROUND"
+echo "========================================================================="
+echo -e "${NC}"
+echo ""
+echo -e "${YELLOW}WARNING:${NC} This script applies a temporary patch to fix a known issue with"
+echo "meson-python that causes duplicate RPATH entries in built Python extensions"
+echo "on macOS, preventing module imports."
+echo ""
+echo "This workaround will:"
+echo "  1. Install fourdst using pip"
+echo "  2. Locate the installed extension module"
+echo "  3. Remove duplicate RPATH entries using install_name_tool"
+echo ""
+echo "This is a temporary solution while the meson-python team resolves the"
+echo "duplicate RPATH bug. For more information, see:"
+echo "  https://github.com/mesonbuild/meson-python/issues/813"
+echo ""
+echo -e "${YELLOW}Do you understand and wish to continue? [y/N]${NC} "
+read -r response
+
+if [[ ! "$response" =~ ^[Yy]$ ]]; then
+    echo -e "${RED}Installation cancelled.${NC}"
+    exit 1
+fi
+
+echo ""
+echo -e "${GREEN}Step 1: Finding current Python environment...${NC}"
+
+# Get the current Python executable
+PYTHON_BIN=$(which python3)
+if [ -z "$PYTHON_BIN" ]; then
+    echo -e "${RED}Error: python3 not found in PATH${NC}"
+    exit 1
+fi
+
+echo "Using Python: $PYTHON_BIN"
+PYTHON_VERSION=$($PYTHON_BIN --version)
+echo "Python version: $PYTHON_VERSION"
+
+# Get site-packages directory
+SITE_PACKAGES=$($PYTHON_BIN -c "import site; print(site.getsitepackages()[0])")
+echo "Site packages: $SITE_PACKAGES"
+echo ""
+
+echo -e "${GREEN}Step 2: Installing fourdst with pip...${NC}"
+$PYTHON_BIN -m pip install . -v
+
+if [ $? -ne 0 ]; then
+    echo -e "${RED}Error: pip install failed${NC}"
+    exit 1
+fi
+echo ""
+
+echo -e "${GREEN}Step 3: Locating installed gridfire extension module...${NC}"
+
+# Find the .so file
+SO_FILE=$(find "$SITE_PACKAGES" -name "gridfire.cpython-*-darwin.so" 2>/dev/null | head -n 1)
+
+if [ -z "$SO_FILE" ]; then
+    echo -e "${RED}Error: Could not find gridfire.cpython-*-darwin.so in $SITE_PACKAGES/gridfire${NC}"
+    echo "Installation may have failed or the file is in an unexpected location."
+    exit 1
+fi
+
+echo "Found gridfire extension module: $SO_FILE"
+echo ""
+
+echo -e "${GREEN}Step 4: Running RPATH fix script for gridfire extension module...${NC}"
+
+# Check if fix_rpath.py exists
+FIX_SCRIPT="build-python/fix_rpaths.py"
+if [ ! -f "$FIX_SCRIPT" ]; then
+    echo -e "${RED}Error: $FIX_SCRIPT not found${NC}"
+    echo "Please ensure you're running this script from the project root directory."
+    exit 1
+fi
+
+# Run the fix script
+$PYTHON_BIN "$FIX_SCRIPT" "$SO_FILE"
+
+if [ $? -ne 0 ]; then
+    echo -e "${RED}Error: RPATH fix script failed for gridfire extension module${NC}"
+    exit 1
+fi
+
+echo -e "${GREEN}Step 5: Locating installed fourdst extension module...${NC}"
+
+# Find the .so file
+SO_FILE=$(find "$SITE_PACKAGES/fourdst" -name "_phys.cpython-*-darwin.so" 2>/dev/null | head -n 1)
+
+if [ -z "$SO_FILE" ]; then
+    echo -e "${RED}Error: Could not find _phys.cpython-*-darwin.so in $SITE_PACKAGES/fourdst${NC}"
+    echo "Installation may have failed or the file is in an unexpected location."
+    exit 1
+fi
+
+echo "Found fourdst extension module: $SO_FILE"
+echo ""
+
+echo -e "${GREEN}Step 6: Running RPATH fix script for fourdst extension module...${NC}"
+
+# Check if fix_rpath.py exists
+FIX_SCRIPT="build-python/fix_rpaths.py"
+if [ ! -f "$FIX_SCRIPT" ]; then
+    echo -e "${RED}Error: $FIX_SCRIPT not found${NC}"
+    echo "Please ensure you're running this script from the project root directory."
+    exit 1
+fi
+
+# Run the fix script
+$PYTHON_BIN "$FIX_SCRIPT" "$SO_FILE"
+
+if [ $? -ne 0 ]; then
+    echo -e "${RED}Error: RPATH fix script failed for fourdst extension module${NC}"
+    exit 1
+fi
+
+echo ""
+echo -e "${GREEN}=========================================================================${NC}"
+echo -e "${GREEN}                    Installation Complete!${NC}"
+echo -e "${GREEN}=========================================================================${NC}"
+echo ""
+echo "You can now use fourdst in your Python environment."
+echo ""
+echo "Test the installation with:"
+echo "  $PYTHON_BIN -c 'import gridfire; print(gridfire.__version__)'"
+echo ""
+echo -e "${YELLOW}Note:${NC} If you reinstall or upgrade fourdst, you will need to run this"
+echo "script again to apply the RPATH fix."
+echo ""
\ No newline at end of file
diff --git a/src/include/gridfire/engine/procedures/construction.h b/src/include/gridfire/engine/procedures/construction.h
index e14c98f5..255df4ea 100644
--- a/src/include/gridfire/engine/procedures/construction.h
+++ b/src/include/gridfire/engine/procedures/construction.h
@@ -22,9 +22,13 @@ namespace gridfire {
         ELECTRON_CAPTURE = 1 << 3, // 8
         POSITRON_CAPTURE = 1 << 4, // 16
 
-        WEAK = BETA_MINUS | BETA_PLUS | ELECTRON_CAPTURE | POSITRON_CAPTURE,
-        DEFAULT = STRONG,
-        ALL = STRONG | WEAK
+        REACLIB_WEAK = 1 << 5,
+
+        WRL_WEAK = BETA_MINUS | BETA_PLUS | ELECTRON_CAPTURE | POSITRON_CAPTURE,
+        REACLIB = STRONG | REACLIB_WEAK,
+        DEFAULT = REACLIB,
+
+        ALL = STRONG | WRL_WEAK
     };
 
     constexpr auto to_underlying(NetworkConstructionFlags f) noexcept {
@@ -45,12 +49,13 @@ namespace gridfire {
 
     inline std::string NetworkConstructionFlagsToString(NetworkConstructionFlags flags) {
         std::stringstream ss;
-        constexpr std::array<NetworkConstructionFlags, 5> bases_flags_array = {
+        constexpr std::array<NetworkConstructionFlags, 6> bases_flags_array = {
             NetworkConstructionFlags::STRONG,
             NetworkConstructionFlags::BETA_MINUS,
             NetworkConstructionFlags::BETA_PLUS,
             NetworkConstructionFlags::ELECTRON_CAPTURE,
-            NetworkConstructionFlags::POSITRON_CAPTURE
+            NetworkConstructionFlags::POSITRON_CAPTURE,
+            NetworkConstructionFlags::REACLIB_WEAK
         };
 
         const std::unordered_map<NetworkConstructionFlags, std::string> bases_string_map = {
@@ -58,7 +63,8 @@ namespace gridfire {
             {NetworkConstructionFlags::BETA_MINUS, "BetaMinus"},
             {NetworkConstructionFlags::BETA_PLUS, "BetaPlus"},
             {NetworkConstructionFlags::ELECTRON_CAPTURE, "ElectronCapture"},
-            {NetworkConstructionFlags::POSITRON_CAPTURE, "PositronCapture"}
+            {NetworkConstructionFlags::POSITRON_CAPTURE, "PositronCapture"},
+            {NetworkConstructionFlags::REACLIB_WEAK, "ReaclibWeak"}
         };
 
         size_t i = 0;
diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp
index 8584e344..faf61de5 100644
--- a/src/lib/engine/procedures/construction.cpp
+++ b/src/lib/engine/procedures/construction.cpp
@@ -44,7 +44,7 @@ namespace {
         const gridfire::NetworkConstructionFlags reactionTypes
     ) {
         gridfire::reaction::ReactionSet weak_reaction_pool;
-        if (!has_flag(reactionTypes, gridfire::NetworkConstructionFlags::WEAK)) {
+        if (!has_flag(reactionTypes, gridfire::NetworkConstructionFlags::WRL_WEAK)) {
             return weak_reaction_pool;
         }
 
@@ -109,13 +109,38 @@ namespace {
         if (has_flag(reaction_types, gridfire::NetworkConstructionFlags::STRONG)) {
             const auto& allReaclibReactions = gridfire::reaclib::get_all_reaclib_reactions();
             for (const auto& reaction : allReaclibReactions) {
-                if (!reaction->is_reverse() && !reaclib_reaction_is_weak(*reaction)) { // Only add reactions of the correct direction and which are not weak. Weak reactions are handled from the WRL separately which provides much higher quality weak reactions than reaclib does
+                const bool isWeakReaction = reaclib_reaction_is_weak(*reaction);
+                const bool okayToUseReaclibWeakReaction = has_flag(reaction_types, gridfire::NetworkConstructionFlags::REACLIB_WEAK);
+
+                const bool reaclibWeakOkay = !isWeakReaction || okayToUseReaclibWeakReaction;
+                if (!reaction->is_reverse() && reaclibWeakOkay) {
                     strong_reaction_pool.add_reaction(reaction->clone());
                 }
             }
         }
         return strong_reaction_pool;
     }
+
+    bool validate_unique_weak_set(gridfire::NetworkConstructionFlags flag) {
+        // This method ensures that weak reactions will only be fetched from either reaclib or the weak reaction library (WRL)
+        // but not both
+        std::array<gridfire::NetworkConstructionFlags, 4> WRL_Flags = {
+            gridfire::NetworkConstructionFlags::BETA_PLUS,
+            gridfire::NetworkConstructionFlags::ELECTRON_CAPTURE,
+            gridfire::NetworkConstructionFlags::POSITRON_CAPTURE,
+            gridfire::NetworkConstructionFlags::BETA_MINUS
+        };
+
+        if (!has_flag(flag, gridfire::NetworkConstructionFlags::REACLIB_WEAK)) {
+            return true;
+        }
+        for (const auto& WRLReactionType : WRL_Flags) {
+            if (has_flag(flag, WRLReactionType)) {
+                return false;
+            }
+        }
+        return true;
+    }
 }
 
 namespace gridfire {
@@ -131,6 +156,11 @@ namespace gridfire {
         NetworkConstructionFlags ReactionTypes
     ) {
         auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        if (!validate_unique_weak_set(ReactionTypes)) {
+            std::string msg = "Cannot construct network since weak reactions from both reaclib and WRL were requested. Only one weak rate source may be used. In network construction this likely means that the flag NetworkConstructionFlags::REACLIB_WEAK was used along with one of the weak flags (BETA_PLUS, BETA_MINUS, ELECTRON_CAPTURE, POSITRON_CAPTURE). These flags are for WRL rates and may not be used in conjunction with reaclib weak rates.";
+            LOG_ERROR(logger, "{}", msg);
+            throw std::logic_error(msg);
+        }
         LOG_INFO(logger, "Constructing network topology from reaction types : {}", NetworkConstructionFlagsToString(ReactionTypes));
 
         if (ReactionTypes == NetworkConstructionFlags::NONE) {
diff --git a/src/python/bindings.cpp b/src/python/bindings.cpp
index 20f89b8f..0d29b21b 100644
--- a/src/python/bindings.cpp
+++ b/src/python/bindings.cpp
@@ -23,30 +23,30 @@ PYBIND11_MODULE(gridfire, m) {
     auto typeMod  = m.def_submodule("type", "GridFire type bindings");
     register_type_bindings(typeMod);
 
-    // auto partitionMod  = m.def_submodule("partition", "GridFire partition function bindings");
-    // register_partition_bindings(partitionMod);
-    //
-    // auto expectationMod  = m.def_submodule("expectations", "GridFire expectations bindings");
-    // register_expectation_bindings(expectationMod);
-    //
-    // auto reactionMod  = m.def_submodule("reaction", "GridFire reaction bindings");
-    // register_reaction_bindings(reactionMod);
-    //
-    // auto screeningMod  = m.def_submodule("screening", "GridFire plasma screening bindings");
-    // register_screening_bindings(screeningMod);
-    //
-    // auto ioMod  = m.def_submodule("io", "GridFire io bindings");
-    // register_io_bindings(ioMod);
-    //
-    // auto exceptionMod  = m.def_submodule("exceptions", "GridFire exceptions bindings");
-    // register_exception_bindings(exceptionMod);
-    //
-    // auto engineMod  = m.def_submodule("engine", "Engine and Engine View bindings");
-    // register_engine_bindings(engineMod);
-    //
-    // auto solverMod  = m.def_submodule("solver", "GridFire numerical solver bindings");
-    // register_solver_bindings(solverMod);
-    //
-    // auto utilsMod  = m.def_submodule("utils", "GridFire utility method bindings");
-    // register_utils_bindings(utilsMod);
+    auto partitionMod  = m.def_submodule("partition", "GridFire partition function bindings");
+    register_partition_bindings(partitionMod);
+
+    auto expectationMod  = m.def_submodule("expectations", "GridFire expectations bindings");
+    register_expectation_bindings(expectationMod);
+
+    auto reactionMod  = m.def_submodule("reaction", "GridFire reaction bindings");
+    register_reaction_bindings(reactionMod);
+
+    auto screeningMod  = m.def_submodule("screening", "GridFire plasma screening bindings");
+    register_screening_bindings(screeningMod);
+
+    auto ioMod  = m.def_submodule("io", "GridFire io bindings");
+    register_io_bindings(ioMod);
+
+    auto exceptionMod  = m.def_submodule("exceptions", "GridFire exceptions bindings");
+    register_exception_bindings(exceptionMod);
+
+    auto engineMod  = m.def_submodule("engine", "Engine and Engine View bindings");
+    register_engine_bindings(engineMod);
+
+    auto solverMod  = m.def_submodule("solver", "GridFire numerical solver bindings");
+    register_solver_bindings(solverMod);
+
+    auto utilsMod  = m.def_submodule("utils", "GridFire utility method bindings");
+    register_utils_bindings(utilsMod);
 }
\ No newline at end of file
diff --git a/src/python/engine/bindings.cpp b/src/python/engine/bindings.cpp
index 7988df0a..c89e8816 100644
--- a/src/python/engine/bindings.cpp
+++ b/src/python/engine/bindings.cpp
@@ -197,11 +197,11 @@ namespace {
 }
 
 void register_engine_bindings(py::module &m) {
+    register_engine_type_bindings(m);
+    register_engine_procedural_bindings(m);
     register_base_engine_bindings(m);
     register_engine_view_bindings(m);
     register_engine_diagnostic_bindings(m);
-    register_engine_procedural_bindings(m);
-    register_engine_type_bindings(m);
 }
 
 void register_base_engine_bindings(const pybind11::module &m) {
diff --git a/src/python/exceptions/bindings.cpp b/src/python/exceptions/bindings.cpp
index 5fe00648..3223fb0c 100644
--- a/src/python/exceptions/bindings.cpp
+++ b/src/python/exceptions/bindings.cpp
@@ -45,14 +45,10 @@ void register_exception_bindings(const py::module &m) {
 
     py::register_exception<gridfire::exceptions::JacobianError>(m, "JacobianError", m.attr("GridFireEngineError"));
 
-    py::register_exception<gridfire::exceptions::StaleJacobianError>(m, "StaleJacobianError", m.attr("JacobianEngineError"));
-    py::register_exception<gridfire::exceptions::UninitializedJacobianError>(m, "UninitializedJacobianError", m.attr("JacobianEngineError"));
-    py::register_exception<gridfire::exceptions::UnknownJacobianError>(m, "UnknownJacobianError", m.attr("JacobianEngineError"));
+    py::register_exception<gridfire::exceptions::StaleJacobianError>(m, "StaleJacobianError", m.attr("JacobianError"));
+    py::register_exception<gridfire::exceptions::UninitializedJacobianError>(m, "UninitializedJacobianError", m.attr("JacobianError"));
+    py::register_exception<gridfire::exceptions::UnknownJacobianError>(m, "UnknownJacobianError", m.attr("JacobianError"));
 
     py::register_exception<gridfire::exceptions::UtilityError>(m, "UtilityError");
     py::register_exception<gridfire::exceptions::HashingError>(m, "HashingError", m.attr("UtilityError"));
-
-
-
-
 }
diff --git a/src/python/meson.build b/src/python/meson.build
index daa460e3..ec84b6be 100644
--- a/src/python/meson.build
+++ b/src/python/meson.build
@@ -1,10 +1,10 @@
 subdir('types')
-#subdir('utils')
-#subdir('expectations')
-#subdir('exceptions')
-#subdir('io')
-#subdir('partition')
-#subdir('reaction')
-#subdir('screening')
-#subdir('engine')
-#subdir('solver')
+subdir('utils')
+subdir('expectations')
+subdir('exceptions')
+subdir('io')
+subdir('partition')
+subdir('reaction')
+subdir('screening')
+subdir('engine')
+subdir('solver')
diff --git a/src/python/solver/bindings.cpp b/src/python/solver/bindings.cpp
index a4035801..3b4fccd9 100644
--- a/src/python/solver/bindings.cpp
+++ b/src/python/solver/bindings.cpp
@@ -1,4 +1,5 @@
 #include <pybind11/pybind11.h>
+#include <pybind11/functional.h> // needed for std::function
 #include <pybind11/stl.h> // Needed for vectors, maps, sets, strings
 #include <pybind11/stl_bind.h> // Needed for binding std::vector, std::map etc. if needed directly
 #include <pybind11/numpy.h>
diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap
index ebf8dfbb..dc7d6db5 100644
--- a/subprojects/fourdst.wrap
+++ b/subprojects/fourdst.wrap
@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/fourdst
-revision = v0.8.4
+revision = v0.8.5
 depth = 1
diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp
index 6bd90bed..16757bdd 100644
--- a/tests/graphnet_sandbox/main.cpp
+++ b/tests/graphnet_sandbox/main.cpp
@@ -115,7 +115,7 @@ int main(int argc, char* argv[]){
     // netIn.tMax = 1e-14;
     netIn.dt0 = 1e-12;
 
-    GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder, NetworkConstructionFlags::STRONG );
+    GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder, NetworkConstructionFlags::DEFAULT);
 
     ReaclibEngine.setPrecomputation(true);
     ReaclibEngine.setUseReverseReactions(false);
diff --git a/validation/pynucastro/GridFireEquiv/evolve.py b/validation/pynucastro/GridFireEquiv/evolve.py
index 6a9640a9..601bd0e1 100644
--- a/validation/pynucastro/GridFireEquiv/evolve.py
+++ b/validation/pynucastro/GridFireEquiv/evolve.py
@@ -1,46 +1,46 @@
-import numpy as np
-from scipy.integrate import solve_ivp
-
-import pp_chain_robust as network
-
-
-T = 1.5e7
-rho = 1.5e2  
-
-Y0 = np.zeros(network.nnuc)
-Y0[network.jp] = 0.702583  
-Y0[network.jhe3] = 9.74903e-6  
-Y0[network.jhe4] = 0.068963
-Y0[network.jc12] = 0.000250029
-Y0[network.jn14] = 7.85632e-5
-Y0[network.jo16] = 0.00060151
-Y0[network.jne20] = 8.10399e-5
-Y0[network.jmg24] = 2.15159e-5
-
-t_start = 0.0
-t_end = 3.14e17
-t_span = (t_start, t_end)
-
-sol = solve_ivp(
-    network.rhs,
-    t_span=t_span,
-    y0=Y0,
-    method='Radau',
-    jac=network.jacobian,
-    args=(rho, T),
-    dense_output=True,
-    rtol=1e-8,
-    atol=1e-8
-)
-
-
-with open("pynucastro_results.csv", 'w') as f:
-    f.write("t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24\n")
-    for (t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24) in zip(sol.t, sol.y[network.jp, :], sol.y[network.jd, :], sol.y[network.jhe3, :], sol.y[network.jhe4, :], sol.y[network.jc12, :], sol.y[network.jn14, :], sol.y[network.jo16, :], sol.y[network.jne20, :], sol.y[network.jmg24, :]):
-        f.write(f"{t},{h1},{h2},{he3},{he4},{c12},{n14},{o16},{ne20},{mg24}\n")
+# import numpy as np
+# from scipy.integrate import solve_ivp
+#
+# import pp_chain_robust as network
+#
+#
+# T = 1e7
+# rho = 1.5e2
+#
+# Y0 = np.zeros(network.nnuc)
+# Y0[network.jp] = 0.702583
+# Y0[network.jhe3] = 9.74903e-6
+# Y0[network.jhe4] = 0.068963
+# Y0[network.jc12] = 0.000250029
+# Y0[network.jn14] = 7.85632e-5
+# Y0[network.jo16] = 0.00060151
+# Y0[network.jne20] = 8.10399e-5
+# Y0[network.jmg24] = 2.15159e-5
+#
+# t_start = 0.0
+# t_end = 3.14e17
+# t_span = (t_start, t_end)
+#
+# sol = solve_ivp(
+#     network.rhs,
+#     t_span=t_span,
+#     y0=Y0,
+#     method='Radau',
+#     jac=network.jacobian,
+#     args=(rho, T),
+#     dense_output=True,
+#     rtol=1e-8,
+#     atol=1e-8
+# )
+#
+#
+# with open("pynucastro_results.csv", 'w') as f:
+#     f.write("t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24\n")
+#     for (t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24) in zip(sol.t, sol.y[network.jp, :], sol.y[network.jd, :], sol.y[network.jhe3, :], sol.y[network.jhe4, :], sol.y[network.jc12, :], sol.y[network.jn14, :], sol.y[network.jo16, :], sol.y[network.jne20, :], sol.y[network.jmg24, :]):
+#         f.write(f"{t},{h1},{h2},{he3},{he4},{c12},{n14},{o16},{ne20},{mg24}\n")
 
 from gridfire.engine import GraphEngine, MultiscalePartitioningEngineView, AdaptiveEngineView
-from gridfire.solver import DirectNetworkSolver
+from gridfire.solver import CVODESolverStrategy
 from gridfire.type import NetIn
 
 from fourdst.composition import Composition
@@ -61,7 +61,7 @@ netIn = NetIn()
 netIn.composition = comp
 netIn.temperature = 1.5e7
 netIn.density = 1.6e2
-netIn.tMax = 3.14e17
+netIn.tMax = 3.14e16
 netIn.dt0 = 1e-12
 
 baseEngine = GraphEngine(netIn.composition, 2)
@@ -71,7 +71,7 @@ qseEngine = MultiscalePartitioningEngineView(baseEngine)
 
 adaptiveEngine = AdaptiveEngineView(qseEngine)
 
-solver = DirectNetworkSolver(adaptiveEngine)
+solver = CVODESolverStrategy(adaptiveEngine)
 
 data = []
 def callback(context):
@@ -98,10 +98,10 @@ def callback(context):
 
 
 solver.set_callback(callback)
-results = solver.evaluate(netIn)
+results = solver.evaluate(netIn, False)
 
-with open("gridfire_results.csv", 'w') as f:
-    f.write("t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24\n")
-    for row in data:
-        rowStr = ','.join([str(x) for x in row])
-        f.write(f"{rowStr}\n")
+# with open("gridfire_results.csv", 'w') as f:
+#     f.write("t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24\n")
+#     for row in data:
+#         rowStr = ','.join([str(x) for x in row])
+#         f.write(f"{rowStr}\n")
diff --git a/validation/pynucastro/GridFireEquiv/generate_proton_proton.py b/validation/pynucastro/GridFireEquiv/generate_proton_proton.py
index 68fbb986..7450f407 100644
--- a/validation/pynucastro/GridFireEquiv/generate_proton_proton.py
+++ b/validation/pynucastro/GridFireEquiv/generate_proton_proton.py
@@ -1,7 +1,6 @@
 import pynucastro as pyna
 
 from gridfire.engine import GraphEngine, MultiscalePartitioningEngineView, AdaptiveEngineView
-from gridfire.solver import DirectNetworkSolver
 from gridfire.type import NetIn
 
 from fourdst.composition import Composition
@@ -19,7 +18,7 @@ print(f"Initial H-1 mass fraction {comp.getMassFraction("H-1")}")
 
 netIn = NetIn()
 netIn.composition = comp
-netIn.temperature = 1.5e7
+netIn.temperature = 1e7
 netIn.density = 1.6e2
 netIn.tMax = 1e-9
 netIn.dt0 = 1e-12
diff --git a/validation/pynucastro/GridFireEquiv/pp_chain_robust.py b/validation/pynucastro/GridFireEquiv/pp_chain_robust.py
index 1482525c..0f817d82 100644
--- a/validation/pynucastro/GridFireEquiv/pp_chain_robust.py
+++ b/validation/pynucastro/GridFireEquiv/pp_chain_robust.py
@@ -12,59 +12,62 @@ jd = 2
 jt = 3
 jhe3 = 4
 jhe4 = 5
-jli6 = 6
-jli7 = 7
-jbe7 = 8
-jb8 = 9
-jc11 = 10
-jc12 = 11
-jc13 = 12
-jn13 = 13
-jn14 = 14
-jn15 = 15
-jo14 = 16
-jo15 = 17
-jo16 = 18
-jo17 = 19
-jo18 = 20
-jf15 = 21
-jf16 = 22
-jf17 = 23
-jf18 = 24
-jf19 = 25
-jne18 = 26
-jne19 = 27
-jne20 = 28
-jne21 = 29
-jna21 = 30
-jna22 = 31
-jna23 = 32
-jna24 = 33
-jmg22 = 34
-jmg23 = 35
-jmg24 = 36
-jmg25 = 37
-jmg26 = 38
-jal25 = 39
-jal26 = 40
-jal27 = 41
-jal28 = 42
-jsi26 = 43
-jsi28 = 44
-jsi29 = 45
-jsi30 = 46
-jp29 = 47
-jp30 = 48
-jp31 = 49
-jp32 = 50
-js31 = 51
-js32 = 52
-js34 = 53
-jcl32 = 54
-jcl34 = 55
-jcl35 = 56
-jar35 = 57
-nnuc = 58
+jhe6 = 6
+jli6 = 7
+jli7 = 8
+jbe7 = 9
+jbe9 = 10
+jb8 = 11
+jc11 = 12
+jc12 = 13
+jc13 = 14
+jc14 = 15
+jn13 = 16
+jn14 = 17
+jn15 = 18
+jo14 = 19
+jo15 = 20
+jo16 = 21
+jo17 = 22
+jo18 = 23
+jf15 = 24
+jf16 = 25
+jf17 = 26
+jf18 = 27
+jf19 = 28
+jne18 = 29
+jne19 = 30
+jne20 = 31
+jne21 = 32
+jna21 = 33
+jna22 = 34
+jna23 = 35
+jna24 = 36
+jmg22 = 37
+jmg23 = 38
+jmg24 = 39
+jmg25 = 40
+jmg26 = 41
+jal25 = 42
+jal26 = 43
+jal27 = 44
+jal28 = 45
+jsi26 = 46
+jsi28 = 47
+jsi29 = 48
+jsi30 = 49
+jp29 = 50
+jp30 = 51
+jp31 = 52
+jp32 = 53
+js31 = 54
+js32 = 55
+js34 = 56
+jcl32 = 57
+jcl34 = 58
+jcl35 = 59
+jar35 = 60
+nnuc = 61
 
 A = np.zeros((nnuc), dtype=np.int32)
 
@@ -74,13 +77,16 @@ A[jd] = 2
 A[jt] = 3
 A[jhe3] = 3
 A[jhe4] = 4
+A[jhe6] = 6
 A[jli6] = 6
 A[jli7] = 7
 A[jbe7] = 7
+A[jbe9] = 9
 A[jb8] = 8
 A[jc11] = 11
 A[jc12] = 12
 A[jc13] = 13
+A[jc14] = 14
 A[jn13] = 13
 A[jn14] = 14
 A[jn15] = 15
@@ -135,13 +141,16 @@ Z[jd] = 1
 Z[jt] = 1
 Z[jhe3] = 2
 Z[jhe4] = 2
+Z[jhe6] = 2
 Z[jli6] = 3
 Z[jli7] = 3
 Z[jbe7] = 4
+Z[jbe9] = 4
 Z[jb8] = 5
 Z[jc11] = 6
 Z[jc12] = 6
 Z[jc13] = 6
+Z[jc14] = 6
 Z[jn13] = 7
 Z[jn14] = 7
 Z[jn15] = 7
@@ -197,13 +206,16 @@ mass[jd] = 0.0030058819195053215
 mass[jt] = 0.004501206494525079
 mass[jhe3] = 0.004501176706825056
 mass[jhe4] = 0.0059735574859708365
+mass[jhe6] = 0.008982694165199753
 mass[jli6] = 0.008977078184259593
 mass[jli7] = 0.010470810414554471
 mass[jbe7] = 0.010472191322584432
+mass[jbe9] = 0.013449944991877256
 mass[jb8] = 0.011976069136782909
 mass[jc11] = 0.016433661161513517
 mass[jc12] = 0.017909017027273523
 mass[jc13] = 0.01940644192976114
+mass[jc14] = 0.020898691600487556
 mass[jn13] = 0.01940999951603316
 mass[jn14] = 0.020898440897976135
 mass[jn15] = 0.022386433805845516
@@ -257,13 +269,16 @@ names.append("H2")
 names.append("H3")
 names.append("He3")
 names.append("He4")
+names.append("He6")
 names.append("Li6")
 names.append("Li7")
 names.append("Be7")
+names.append("Be9")
 names.append("B8")
 names.append("C11")
 names.append("C12")
 names.append("C13")
+names.append("C14")
 names.append("N13")
 names.append("N14")
 names.append("N15")
@@ -332,7 +347,9 @@ def energy_release(dY):
     ("n__p__weak__wc12", numba.float64),
     ("t__He3__weak__wc12", numba.float64),
     ("He3__t__weak__electron_capture", numba.float64),
+    ("He6__Li6__weak__wc12", numba.float64),
     ("Be7__Li7__weak__electron_capture", numba.float64),
+    ("C14__N14__weak__wc12", numba.float64),
     ("N13__C13__weak__wc12", numba.float64),
     ("O14__N14__weak__wc12", numba.float64),
     ("O15__N15__weak__wc12", numba.float64),
@@ -371,10 +388,12 @@ def energy_release(dY):
     ("C11__He4_Be7", numba.float64),
     ("C12__n_C11", numba.float64),
     ("C13__n_C12", numba.float64),
+    ("C14__n_C13", numba.float64),
     ("N13__p_C12", numba.float64),
     ("N14__n_N13", numba.float64),
     ("N14__p_C13", numba.float64),
     ("N15__n_N14", numba.float64),
+    ("N15__p_C14", numba.float64),
     ("O14__p_N13", numba.float64),
     ("O15__n_O14", numba.float64),
     ("O15__p_N14", numba.float64),
@@ -383,6 +402,7 @@ def energy_release(dY):
     ("O16__He4_C12", numba.float64),
     ("O17__n_O16", numba.float64),
     ("O18__n_O17", numba.float64),
+    ("O18__He4_C14", numba.float64),
     ("F15__p_O14", numba.float64),
     ("F16__n_F15", numba.float64),
     ("F16__p_O15", numba.float64),
@@ -467,7 +487,9 @@ def energy_release(dY):
     ("Cl35__He4_P31", numba.float64),
     ("Ar35__p_Cl34", numba.float64),
     ("Ar35__He4_S31", numba.float64),
+    ("He6__n_n_He4", numba.float64),
     ("Li6__n_p_He4", numba.float64),
+    ("Be9__n_He4_He4", numba.float64),
     ("C12__He4_He4_He4", numba.float64),
     ("n_p__d", numba.float64),
     ("p_p__d__weak__bet_pos_", numba.float64),
@@ -489,7 +511,10 @@ def energy_release(dY):
     ("n_C12__C13", numba.float64),
     ("p_C12__N13", numba.float64),
     ("He4_C12__O16", numba.float64),
+    ("n_C13__C14", numba.float64),
     ("p_C13__N14", numba.float64),
+    ("p_C14__N15", numba.float64),
+    ("He4_C14__O18", numba.float64),
     ("n_N13__N14", numba.float64),
     ("p_N13__O14", numba.float64),
     ("n_N14__N15", numba.float64),
@@ -601,16 +626,22 @@ def energy_release(dY):
     ("p_Li6__He4_He3", numba.float64),
     ("d_Li6__n_Be7", numba.float64),
     ("d_Li6__p_Li7", numba.float64),
+    ("He4_Li6__p_Be9", numba.float64),
     ("p_Li7__n_Be7", numba.float64),
     ("p_Li7__d_Li6", numba.float64),
     ("p_Li7__He4_He4", numba.float64),
+    ("t_Li7__n_Be9", numba.float64),
     ("n_Be7__p_Li7", numba.float64),
     ("n_Be7__d_Li6", numba.float64),
     ("n_Be7__He4_He4", numba.float64),
+    ("n_Be9__t_Li7", numba.float64),
+    ("p_Be9__He4_Li6", numba.float64),
+    ("He4_Be9__n_C12", numba.float64),
     ("He4_B8__p_C11", numba.float64),
     ("p_C11__He4_B8", numba.float64),
     ("He4_C11__n_O14", numba.float64),
     ("He4_C11__p_N14", numba.float64),
+    ("n_C12__He4_Be9", numba.float64),
     ("He4_C12__n_O15", numba.float64),
     ("He4_C12__p_N15", numba.float64),
     ("C12_C12__n_Mg23", numba.float64),
@@ -619,14 +650,19 @@ def energy_release(dY):
     ("p_C13__n_N13", numba.float64),
     ("d_C13__n_N14", numba.float64),
     ("He4_C13__n_O16", numba.float64),
+    ("p_C14__n_N14", numba.float64),
+    ("d_C14__n_N15", numba.float64),
+    ("He4_C14__n_O17", numba.float64),
     ("n_N13__p_C13", numba.float64),
     ("He4_N13__n_F16", numba.float64),
     ("He4_N13__p_O16", numba.float64),
+    ("n_N14__p_C14", numba.float64),
     ("n_N14__d_C13", numba.float64),
     ("p_N14__n_O14", numba.float64),
     ("p_N14__He4_C11", numba.float64),
     ("He4_N14__n_F17", numba.float64),
     ("He4_N14__p_O17", numba.float64),
+    ("n_N15__d_C14", numba.float64),
     ("p_N15__n_O15", numba.float64),
     ("p_N15__He4_C12", numba.float64),
     ("He4_N15__n_F18", numba.float64),
@@ -649,6 +685,7 @@ def energy_release(dY):
     ("O16_O16__n_S31", numba.float64),
     ("O16_O16__p_P31", numba.float64),
     ("O16_O16__He4_Si28", numba.float64),
+    ("n_O17__He4_C14", numba.float64),
     ("p_O17__n_F17", numba.float64),
     ("p_O17__He4_N14", numba.float64),
     ("He4_O17__n_Ne20", numba.float64),
@@ -815,13 +852,17 @@ def energy_release(dY):
     ("He3_He3__p_p_He4", numba.float64),
     ("d_Li7__n_He4_He4", numba.float64),
     ("d_Be7__p_He4_He4", numba.float64),
+    ("p_Be9__d_He4_He4", numba.float64),
     ("n_B8__p_He4_He4", numba.float64),
     ("n_C11__He4_He4_He4", numba.float64),
     ("t_Li7__n_n_He4_He4", numba.float64),
     ("He3_Li7__n_p_He4_He4", numba.float64),
     ("t_Be7__n_p_He4_He4", numba.float64),
     ("He3_Be7__p_p_He4_He4", numba.float64),
+    ("p_Be9__n_p_He4_He4", numba.float64),
+    ("n_n_He4__He6", numba.float64),
     ("n_p_He4__Li6", numba.float64),
+    ("n_He4_He4__Be9", numba.float64),
     ("He4_He4_He4__C12", numba.float64),
     ("n_p_p__p_d", numba.float64),
     ("n_n_He4__t_t", numba.float64),
@@ -830,10 +871,12 @@ def energy_release(dY):
     ("n_He4_He4__d_Li7", numba.float64),
     ("p_He4_He4__n_B8", numba.float64),
     ("p_He4_He4__d_Be7", numba.float64),
+    ("d_He4_He4__p_Be9", numba.float64),
     ("He4_He4_He4__n_C11", numba.float64),
     ("n_n_He4_He4__t_Li7", numba.float64),
     ("n_p_He4_He4__He3_Li7", numba.float64),
     ("n_p_He4_He4__t_Be7", numba.float64),
+    ("n_p_He4_He4__p_Be9", numba.float64),
     ("p_p_He4_He4__He3_Be7", numba.float64),
 ])
 class RateEval:
@@ -841,7 +884,9 @@ class RateEval:
         self.n__p__weak__wc12 = np.nan
         self.t__He3__weak__wc12 = np.nan
         self.He3__t__weak__electron_capture = np.nan
+        self.He6__Li6__weak__wc12 = np.nan
         self.Be7__Li7__weak__electron_capture = np.nan
+        self.C14__N14__weak__wc12 = np.nan
         self.N13__C13__weak__wc12 = np.nan
         self.O14__N14__weak__wc12 = np.nan
         self.O15__N15__weak__wc12 = np.nan
@@ -880,10 +925,12 @@ class RateEval:
         self.C11__He4_Be7 = np.nan
         self.C12__n_C11 = np.nan
         self.C13__n_C12 = np.nan
+        self.C14__n_C13 = np.nan
         self.N13__p_C12 = np.nan
         self.N14__n_N13 = np.nan
         self.N14__p_C13 = np.nan
         self.N15__n_N14 = np.nan
+        self.N15__p_C14 = np.nan
         self.O14__p_N13 = np.nan
         self.O15__n_O14 = np.nan
         self.O15__p_N14 = np.nan
@@ -892,6 +939,7 @@ class RateEval:
         self.O16__He4_C12 = np.nan
         self.O17__n_O16 = np.nan
         self.O18__n_O17 = np.nan
+        self.O18__He4_C14 = np.nan
         self.F15__p_O14 = np.nan
         self.F16__n_F15 = np.nan
         self.F16__p_O15 = np.nan
@@ -976,7 +1024,9 @@ class RateEval:
         self.Cl35__He4_P31 = np.nan
         self.Ar35__p_Cl34 = np.nan
         self.Ar35__He4_S31 = np.nan
+        self.He6__n_n_He4 = np.nan
         self.Li6__n_p_He4 = np.nan
+        self.Be9__n_He4_He4 = np.nan
         self.C12__He4_He4_He4 = np.nan
         self.n_p__d = np.nan
         self.p_p__d__weak__bet_pos_ = np.nan
@@ -998,7 +1048,10 @@ class RateEval:
         self.n_C12__C13 = np.nan
         self.p_C12__N13 = np.nan
         self.He4_C12__O16 = np.nan
+        self.n_C13__C14 = np.nan
         self.p_C13__N14 = np.nan
+        self.p_C14__N15 = np.nan
+        self.He4_C14__O18 = np.nan
         self.n_N13__N14 = np.nan
         self.p_N13__O14 = np.nan
         self.n_N14__N15 = np.nan
@@ -1110,16 +1163,22 @@ class RateEval:
         self.p_Li6__He4_He3 = np.nan
         self.d_Li6__n_Be7 = np.nan
         self.d_Li6__p_Li7 = np.nan
+        self.He4_Li6__p_Be9 = np.nan
         self.p_Li7__n_Be7 = np.nan
         self.p_Li7__d_Li6 = np.nan
         self.p_Li7__He4_He4 = np.nan
+        self.t_Li7__n_Be9 = np.nan
         self.n_Be7__p_Li7 = np.nan
         self.n_Be7__d_Li6 = np.nan
         self.n_Be7__He4_He4 = np.nan
+        self.n_Be9__t_Li7 = np.nan
+        self.p_Be9__He4_Li6 = np.nan
+        self.He4_Be9__n_C12 = np.nan
         self.He4_B8__p_C11 = np.nan
         self.p_C11__He4_B8 = np.nan
         self.He4_C11__n_O14 = np.nan
         self.He4_C11__p_N14 = np.nan
+        self.n_C12__He4_Be9 = np.nan
         self.He4_C12__n_O15 = np.nan
         self.He4_C12__p_N15 = np.nan
         self.C12_C12__n_Mg23 = np.nan
@@ -1128,14 +1187,19 @@ class RateEval:
         self.p_C13__n_N13 = np.nan
         self.d_C13__n_N14 = np.nan
         self.He4_C13__n_O16 = np.nan
+        self.p_C14__n_N14 = np.nan
+        self.d_C14__n_N15 = np.nan
+        self.He4_C14__n_O17 = np.nan
         self.n_N13__p_C13 = np.nan
         self.He4_N13__n_F16 = np.nan
         self.He4_N13__p_O16 = np.nan
+        self.n_N14__p_C14 = np.nan
         self.n_N14__d_C13 = np.nan
         self.p_N14__n_O14 = np.nan
         self.p_N14__He4_C11 = np.nan
         self.He4_N14__n_F17 = np.nan
         self.He4_N14__p_O17 = np.nan
+        self.n_N15__d_C14 = np.nan
         self.p_N15__n_O15 = np.nan
         self.p_N15__He4_C12 = np.nan
         self.He4_N15__n_F18 = np.nan
@@ -1158,6 +1222,7 @@ class RateEval:
         self.O16_O16__n_S31 = np.nan
         self.O16_O16__p_P31 = np.nan
         self.O16_O16__He4_Si28 = np.nan
+        self.n_O17__He4_C14 = np.nan
         self.p_O17__n_F17 = np.nan
         self.p_O17__He4_N14 = np.nan
         self.He4_O17__n_Ne20 = np.nan
@@ -1324,13 +1389,17 @@ class RateEval:
         self.He3_He3__p_p_He4 = np.nan
         self.d_Li7__n_He4_He4 = np.nan
         self.d_Be7__p_He4_He4 = np.nan
+        self.p_Be9__d_He4_He4 = np.nan
         self.n_B8__p_He4_He4 = np.nan
         self.n_C11__He4_He4_He4 = np.nan
         self.t_Li7__n_n_He4_He4 = np.nan
         self.He3_Li7__n_p_He4_He4 = np.nan
         self.t_Be7__n_p_He4_He4 = np.nan
         self.He3_Be7__p_p_He4_He4 = np.nan
+        self.p_Be9__n_p_He4_He4 = np.nan
+        self.n_n_He4__He6 = np.nan
         self.n_p_He4__Li6 = np.nan
+        self.n_He4_He4__Be9 = np.nan
         self.He4_He4_He4__C12 = np.nan
         self.n_p_p__p_d = np.nan
         self.n_n_He4__t_t = np.nan
@@ -1339,10 +1408,12 @@ class RateEval:
         self.n_He4_He4__d_Li7 = np.nan
         self.p_He4_He4__n_B8 = np.nan
         self.p_He4_He4__d_Be7 = np.nan
+        self.d_He4_He4__p_Be9 = np.nan
         self.He4_He4_He4__n_C11 = np.nan
         self.n_n_He4_He4__t_Li7 = np.nan
         self.n_p_He4_He4__He3_Li7 = np.nan
         self.n_p_He4_He4__t_Be7 = np.nan
+        self.n_p_He4_He4__p_Be9 = np.nan
         self.p_p_He4_He4__He3_Be7 = np.nan
 
 @numba.njit()
@@ -1380,6 +1451,16 @@ def He3__t__weak__electron_capture(rate_eval, tf):
 
     rate_eval.He3__t__weak__electron_capture = rate
 
+@numba.njit()
+def He6__Li6__weak__wc12(rate_eval, tf):
+    # He6 --> Li6
+    rate = 0.0
+
+    # wc12w
+    rate += np.exp(  -0.144619)
+
+    rate_eval.He6__Li6__weak__wc12 = rate
+
 @numba.njit()
 def Be7__Li7__weak__electron_capture(rate_eval, tf):
     # Be7 --> Li7
@@ -1391,6 +1472,16 @@ def Be7__Li7__weak__electron_capture(rate_eval, tf):
 
     rate_eval.Be7__Li7__weak__electron_capture = rate
 
+@numba.njit()
+def C14__N14__weak__wc12(rate_eval, tf):
+    # C14 --> N14
+    rate = 0.0
+
+    # wc12w
+    rate += np.exp(  -26.2827)
+
+    rate_eval.C14__N14__weak__wc12 = rate
+
 @numba.njit()
 def N13__C13__weak__wc12(rate_eval, tf):
     # N13 --> C13
@@ -1820,6 +1911,17 @@ def C13__n_C12(rate_eval, tf):
 
     rate_eval.C13__n_C12 = rate
 
+@numba.njit()
+def C14__n_C13(rate_eval, tf):
+    # C14 --> n + C13
+    rate = 0.0
+
+    # ks03 
+    rate += np.exp(  59.5926 + -95.0156*tf.T9i + 18.3578*tf.T913i + -46.5786*tf.T913
+                  + 2.58472*tf.T9 + -0.118622*tf.T953 + 21.4142*tf.lnT9)
+
+    rate_eval.C14__n_C13 = rate
+
 @numba.njit()
 def N13__p_C12(rate_eval, tf):
     # N13 --> p + C12
@@ -1875,6 +1977,20 @@ def N15__n_N14(rate_eval, tf):
 
     rate_eval.N15__n_N14 = rate
 
+@numba.njit()
+def N15__p_C14(rate_eval, tf):
+    # N15 --> p + C14
+    rate = 0.0
+
+    # il10r
+    rate += np.exp(  40.0115 + -119.975*tf.T9i + -10.658*tf.T913i + 1.73644*tf.T913
+                  + -0.350498*tf.T9 + 0.0279902*tf.T953)
+    # il10n
+    rate += np.exp(  43.0281 + -118.452*tf.T9i + -13.9619*tf.T913i + -4.34315*tf.T913
+                  + 6.64922*tf.T9 + -3.22592*tf.T953 + 0.833333*tf.lnT9)
+
+    rate_eval.N15__p_C14 = rate
+
 @numba.njit()
 def O14__p_N13(rate_eval, tf):
     # O14 --> p + N13
@@ -1984,6 +2100,22 @@ def O18__n_O17(rate_eval, tf):
 
     rate_eval.O18__n_O17 = rate
 
+@numba.njit()
+def O18__He4_C14(rate_eval, tf):
+    # O18 --> He4 + C14
+    rate = 0.0
+
+    # il10r
+    rate += np.exp(  36.546 + -82.6514*tf.T9i + -3.83188*tf.T913
+                  + 1.64358*tf.T9 + -0.177785*tf.T953)
+    # il10r
+    rate += np.exp(  0.910093 + -74.3219*tf.T9i)
+    # il10n
+    rate += np.exp(  43.2028 + -72.2531*tf.T9i + -31.7222*tf.T913i + 11.3923*tf.T913
+                  + -9.92249*tf.T9 + -2.0*tf.T953 + 0.833333*tf.lnT9)
+
+    rate_eval.O18__He4_C14 = rate
+
 @numba.njit()
 def F15__p_O14(rate_eval, tf):
     # F15 --> p + O14
@@ -3103,6 +3235,17 @@ def Ar35__He4_S31(rate_eval, tf):
 
     rate_eval.Ar35__He4_S31 = rate
 
+@numba.njit()
+def He6__n_n_He4(rate_eval, tf):
+    # He6 --> n + n + He4
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  22.178 + -20.8994*tf.T9i + 0.694279*tf.T913i + -3.33326*tf.T913
+                  + 0.507932*tf.T9 + -0.0427342*tf.T953 + 2.0*tf.lnT9)
+
+    rate_eval.He6__n_n_He4 = rate
+
 @numba.njit()
 def Li6__n_p_He4(rate_eval, tf):
     # Li6 --> n + p + He4
@@ -3114,6 +3257,20 @@ def Li6__n_p_He4(rate_eval, tf):
 
     rate_eval.Li6__n_p_He4 = rate
 
+@numba.njit()
+def Be9__n_He4_He4(rate_eval, tf):
+    # Be9 --> n + He4 + He4
+    rate = 0.0
+
+    # ac12n
+    rate += np.exp(  37.273 + -18.2597*tf.T9i + -13.3317*tf.T913i + 13.2237*tf.T913
+                  + -9.06339*tf.T9 + 2.33333*tf.lnT9)
+    # ac12r
+    rate += np.exp(  38.6902 + -19.2792*tf.T9i + -1.56673*tf.T913i + -5.43497*tf.T913
+                  + 0.673807*tf.T9 + -0.041014*tf.T953 + 1.5*tf.lnT9)
+
+    rate_eval.Be9__n_He4_He4 = rate
+
 @numba.njit()
 def C12__He4_He4_He4(rate_eval, tf):
     # C12 --> 3 He4
@@ -3391,6 +3548,17 @@ def He4_C12__O16(rate_eval, tf):
 
     rate_eval.He4_C12__O16 = rate
 
+@numba.njit()
+def n_C13__C14(rate_eval, tf):
+    # C13 + n --> C14
+    rate = 0.0
+
+    # ks03 
+    rate += np.exp(  35.3048 + -0.133687*tf.T9i + 18.3578*tf.T913i + -46.5786*tf.T913
+                  + 2.58472*tf.T9 + -0.118622*tf.T953 + 19.9142*tf.lnT9)
+
+    rate_eval.n_C13__C14 = rate
+
 @numba.njit()
 def p_C13__N14(rate_eval, tf):
     # C13 + p --> N14
@@ -3408,6 +3576,37 @@ def p_C13__N14(rate_eval, tf):
 
     rate_eval.p_C13__N14 = rate
 
+@numba.njit()
+def p_C14__N15(rate_eval, tf):
+    # C14 + p --> N15
+    rate = 0.0
+
+    # il10r
+    rate += np.exp(  17.1024 + -1.52341*tf.T9i + -10.658*tf.T913i + 1.73644*tf.T913
+                  + -0.350498*tf.T9 + 0.0279902*tf.T953 + -1.5*tf.lnT9)
+    # il10n
+    rate += np.exp(  20.119 + -13.9619*tf.T913i + -4.34315*tf.T913
+                  + 6.64922*tf.T9 + -3.22592*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.p_C14__N15 = rate
+
+@numba.njit()
+def He4_C14__O18(rate_eval, tf):
+    # C14 + He4 --> O18
+    rate = 0.0
+
+    # il10r
+    rate += np.exp(  11.8309 + -10.3983*tf.T9i + -3.83188*tf.T913
+                  + 1.64358*tf.T9 + -0.177785*tf.T953 + -1.5*tf.lnT9)
+    # il10r
+    rate += np.exp(  -23.805 + -2.06876*tf.T9i
+                  + -1.5*tf.lnT9)
+    # il10n
+    rate += np.exp(  18.4877 + -31.7222*tf.T913i + 11.3923*tf.T913
+                  + -9.92249*tf.T9 + -2.0*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.He4_C14__O18 = rate
+
 @numba.njit()
 def n_N13__N14(rate_eval, tf):
     # N13 + n --> N14
@@ -4912,6 +5111,23 @@ def d_Li6__p_Li7(rate_eval, tf):
 
     rate_eval.d_Li6__p_Li7 = rate
 
+@numba.njit()
+def He4_Li6__p_Be9(rate_eval, tf):
+    # Li6 + He4 --> p + Be9
+    rate = 0.0
+
+    # cf88n
+    rate += np.exp(  25.5847 + -24.6712*tf.T9i + -10.359*tf.T913i + 0.102577*tf.T913
+                  + 4.43544*tf.T9 + -5.97105*tf.T953 + -0.666667*tf.lnT9)
+    # cf88r
+    rate += np.exp(  19.8324 + -29.8312*tf.T9i
+                  + -0.75*tf.lnT9)
+    # cf88r
+    rate += np.exp(  19.4366 + -27.7172*tf.T9i
+                  + -1.0*tf.lnT9)
+
+    rate_eval.He4_Li6__p_Be9 = rate
+
 @numba.njit()
 def p_Li7__n_Be7(rate_eval, tf):
     # Li7 + p --> n + Be7
@@ -4954,6 +5170,20 @@ def p_Li7__He4_He4(rate_eval, tf):
 
     rate_eval.p_Li7__He4_He4 = rate
 
+@numba.njit()
+def t_Li7__n_Be9(rate_eval, tf):
+    # Li7 + t --> n + Be9
+    rate = 0.0
+
+    # bk91 
+    rate += np.exp(  30.2619 + -11.333*tf.T913i + -0.170964*tf.T913
+                  + -6.30572*tf.T9 + 1.2248*tf.T953 + -0.666667*tf.lnT9)
+    # bk91 
+    rate += np.exp(  34.634 + -11.333*tf.T913i + -7.3964*tf.T913
+                  + 0.947759*tf.T9 + -0.0839008*tf.T953 + 0.333333*tf.lnT9)
+
+    rate_eval.t_Li7__n_Be9 = rate
+
 @numba.njit()
 def n_Be7__p_Li7(rate_eval, tf):
     # Be7 + n --> p + Li7
@@ -4987,6 +5217,59 @@ def n_Be7__He4_He4(rate_eval, tf):
 
     rate_eval.n_Be7__He4_He4 = rate
 
+@numba.njit()
+def n_Be9__t_Li7(rate_eval, tf):
+    # Be9 + n --> t + Li7
+    rate = 0.0
+
+    # bk91 
+    rate += np.exp(  35.9049 + -121.123*tf.T9i + -11.333*tf.T913i + -7.3964*tf.T913
+                  + 0.947759*tf.T9 + -0.0839008*tf.T953 + 0.333333*tf.lnT9)
+    # bk91 
+    rate += np.exp(  31.5328 + -121.123*tf.T9i + -11.333*tf.T913i + -0.170964*tf.T913
+                  + -6.30572*tf.T9 + 1.2248*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.n_Be9__t_Li7 = rate
+
+@numba.njit()
+def p_Be9__He4_Li6(rate_eval, tf):
+    # Be9 + p --> He4 + Li6
+    rate = 0.0
+
+    # cf88n
+    rate += np.exp(  26.0751 + -10.359*tf.T913i + 0.102577*tf.T913
+                  + 4.43544*tf.T9 + -5.97105*tf.T953 + -0.666667*tf.lnT9)
+    # cf88r
+    rate += np.exp(  20.3228 + -5.16*tf.T9i
+                  + -0.75*tf.lnT9)
+    # cf88r
+    rate += np.exp(  19.927 + -3.046*tf.T9i
+                  + -1.0*tf.lnT9)
+
+    rate_eval.p_Be9__He4_Li6 = rate
+
+@numba.njit()
+def He4_Be9__n_C12(rate_eval, tf):
+    # Be9 + He4 --> n + C12
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  11.744 + -4.179*tf.T9i
+                  + -1.5*tf.lnT9)
+    # cf88r
+    rate += np.exp(  -1.48281 + -1.834*tf.T9i
+                  + -1.5*tf.lnT9)
+    # cf88r
+    rate += np.exp(  -9.51959 + -1.184*tf.T9i
+                  + -1.5*tf.lnT9)
+    # cf88n
+    rate += np.exp(  31.464 + -23.87*tf.T913i + 0.566698*tf.T913
+                  + 44.0957*tf.T9 + -314.232*tf.T953 + -0.666667*tf.lnT9)
+    # cf88r
+    rate += np.exp(  19.2962 + -12.732*tf.T9i)
+
+    rate_eval.He4_Be9__n_C12 = rate
+
 @numba.njit()
 def He4_B8__p_C11(rate_eval, tf):
     # B8 + He4 --> p + C11
@@ -5037,6 +5320,28 @@ def He4_C11__p_N14(rate_eval, tf):
 
     rate_eval.He4_C11__p_N14 = rate
 
+@numba.njit()
+def n_C12__He4_Be9(rate_eval, tf):
+    # C12 + n --> He4 + Be9
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  0.858256 + -67.9913*tf.T9i
+                  + -1.5*tf.lnT9)
+    # cf88r
+    rate += np.exp(  -7.17852 + -67.3413*tf.T9i
+                  + -1.5*tf.lnT9)
+    # cf88n
+    rate += np.exp(  33.8051 + -66.1573*tf.T9i + -23.87*tf.T913i + 0.566698*tf.T913
+                  + 44.0957*tf.T9 + -314.232*tf.T953 + -0.666667*tf.lnT9)
+    # cf88r
+    rate += np.exp(  21.6373 + -78.8893*tf.T9i)
+    # cf88r
+    rate += np.exp(  14.0851 + -70.3363*tf.T9i
+                  + -1.5*tf.lnT9)
+
+    rate_eval.n_C12__He4_Be9 = rate
+
 @numba.njit()
 def He4_C12__n_O15(rate_eval, tf):
     # C12 + He4 --> n + O15
@@ -5137,6 +5442,50 @@ def He4_C13__n_O16(rate_eval, tf):
 
     rate_eval.He4_C13__n_O16 = rate
 
+@numba.njit()
+def p_C14__n_N14(rate_eval, tf):
+    # C14 + p --> n + N14
+    rate = 0.0
+
+    # cf88 
+    rate += np.exp(  5.23589 + -7.26442*tf.T9i + 12.3428*tf.T913
+                  + -2.70025*tf.T9 + 0.236625*tf.T953 + 1.0*tf.lnT9)
+    # cf88n
+    rate += np.exp(  14.7608 + -7.26442*tf.T9i + -4.33989*tf.T913
+                  + 11.4311*tf.T9 + -11.7764*tf.T953)
+
+    rate_eval.p_C14__n_N14 = rate
+
+@numba.njit()
+def d_C14__n_N15(rate_eval, tf):
+    # C14 + d --> n + N15
+    rate = 0.0
+
+    # bk92 
+    rate += np.exp(  33.5637 + -16.939*tf.T913i + -4.14392*tf.T913
+                  + 0.438623*tf.T9 + -0.0354193*tf.T953 + 0.333333*tf.lnT9)
+    # bk92 
+    rate += np.exp(  30.6841 + -16.939*tf.T913i + -0.582342*tf.T913
+                  + -8.17066*tf.T9 + 1.70865*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.d_C14__n_N15 = rate
+
+@numba.njit()
+def He4_C14__n_O17(rate_eval, tf):
+    # C14 + He4 --> n + O17
+    rate = 0.0
+
+    # kg91r
+    rate += np.exp(  18.7889 + -23.8496*tf.T9i
+                  + -1.5*tf.lnT9)
+    # kg91r
+    rate += np.exp(  14.4267 + -23.0516*tf.T9i
+                  + -1.5*tf.lnT9)
+    # kg91n
+    rate += np.exp(  11.0417 + -21.0906*tf.T9i)
+
+    rate_eval.He4_C14__n_O17 = rate
+
 @numba.njit()
 def n_N13__p_C13(rate_eval, tf):
     # N13 + n --> p + C13
@@ -5170,6 +5519,20 @@ def He4_N13__p_O16(rate_eval, tf):
 
     rate_eval.He4_N13__p_O16 = rate
 
+@numba.njit()
+def n_N14__p_C14(rate_eval, tf):
+    # N14 + n --> p + C14
+    rate = 0.0
+
+    # cf88 
+    rate += np.exp(  4.13728 + 12.3428*tf.T913
+                  + -2.70025*tf.T9 + 0.236625*tf.T953 + 1.0*tf.lnT9)
+    # cf88n
+    rate += np.exp(  13.6622 + -4.33989*tf.T913
+                  + 11.4311*tf.T9 + -11.7764*tf.T953)
+
+    rate_eval.n_N14__p_C14 = rate
+
 @numba.njit()
 def n_N14__d_C13(rate_eval, tf):
     # N14 + n --> d + C13
@@ -5234,6 +5597,20 @@ def He4_N14__p_O17(rate_eval, tf):
 
     rate_eval.He4_N14__p_O17 = rate
 
+@numba.njit()
+def n_N15__d_C14(rate_eval, tf):
+    # N15 + n --> d + C14
+    rate = 0.0
+
+    # bk92 
+    rate += np.exp(  34.2122 + -92.6344*tf.T9i + -16.939*tf.T913i + -4.14392*tf.T913
+                  + 0.438623*tf.T9 + -0.0354193*tf.T953 + 0.333333*tf.lnT9)
+    # bk92 
+    rate += np.exp(  31.3326 + -92.6344*tf.T9i + -16.939*tf.T913i + -0.582342*tf.T913
+                  + -8.17066*tf.T9 + 1.70865*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.n_N15__d_C14 = rate
+
 @numba.njit()
 def p_N15__n_O15(rate_eval, tf):
     # N15 + p --> n + O15
@@ -5530,6 +5907,22 @@ def O16_O16__He4_Si28(rate_eval, tf):
 
     rate_eval.O16_O16__He4_Si28 = rate
 
+@numba.njit()
+def n_O17__He4_C14(rate_eval, tf):
+    # O17 + n --> He4 + C14
+    rate = 0.0
+
+    # kg91r
+    rate += np.exp(  13.73 + -1.961*tf.T9i
+                  + -1.5*tf.lnT9)
+    # kg91n
+    rate += np.exp(  10.345)
+    # kg91r
+    rate += np.exp(  18.0922 + -2.759*tf.T9i
+                  + -1.5*tf.lnT9)
+
+    rate_eval.n_O17__He4_C14 = rate
+
 @numba.njit()
 def p_O17__n_F17(rate_eval, tf):
     # O17 + p --> n + F17
@@ -7560,6 +7953,23 @@ def d_Be7__p_He4_He4(rate_eval, tf):
 
     rate_eval.d_Be7__p_He4_He4 = rate
 
+@numba.njit()
+def p_Be9__d_He4_He4(rate_eval, tf):
+    # Be9 + p --> d + He4 + He4
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  20.5607 + -5.8*tf.T9i
+                  + -0.75*tf.lnT9)
+    # cf88r
+    rate += np.exp(  20.1768 + -3.046*tf.T9i
+                  + -1.0*tf.lnT9)
+    # cf88n
+    rate += np.exp(  26.0751 + -10.359*tf.T913i + 0.103955*tf.T913
+                  + 4.4262*tf.T9 + -5.95664*tf.T953 + -0.666667*tf.lnT9)
+
+    rate_eval.p_Be9__d_He4_He4 = rate
+
 @numba.njit()
 def n_B8__p_He4_He4(rate_eval, tf):
     # B8 + n --> p + He4 + He4
@@ -7631,6 +8041,31 @@ def He3_Be7__p_p_He4_He4(rate_eval, tf):
 
     rate_eval.He3_Be7__p_p_He4_He4 = rate
 
+@numba.njit()
+def p_Be9__n_p_He4_He4(rate_eval, tf):
+    # Be9 + p --> n + p + He4 + He4
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  20.7431 + -26.725*tf.T9i + 1.40505e-06*tf.T913i + -1.47128e-06*tf.T913
+                  + 6.89313e-08*tf.T9 + -3.55179e-09*tf.T953 + -1.5*tf.lnT9)
+    # cf88n
+    rate += np.exp(  14.6035 + -21.4742*tf.T9i + -0.634849*tf.T913i + 4.82033*tf.T913
+                  + -0.257317*tf.T9 + 0.0134206*tf.T953 + -1.08885*tf.lnT9)
+
+    rate_eval.p_Be9__n_p_He4_He4 = rate
+
+@numba.njit()
+def n_n_He4__He6(rate_eval, tf):
+    # n + n + He4 --> He6
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  -23.9322 + -9.585*tf.T9i + 0.694279*tf.T913i + -3.33326*tf.T913
+                  + 0.507932*tf.T9 + -0.0427342*tf.T953 + -1.0*tf.lnT9)
+
+    rate_eval.n_n_He4__He6 = rate
+
 @numba.njit()
 def n_p_He4__Li6(rate_eval, tf):
     # n + p + He4 --> Li6
@@ -7642,6 +8077,20 @@ def n_p_He4__Li6(rate_eval, tf):
 
     rate_eval.n_p_He4__Li6 = rate
 
+@numba.njit()
+def n_He4_He4__Be9(rate_eval, tf):
+    # n + He4 + He4 --> Be9
+    rate = 0.0
+
+    # ac12r
+    rate += np.exp(  -6.81178 + -1.01953*tf.T9i + -1.56673*tf.T913i + -5.43497*tf.T913
+                  + 0.673807*tf.T9 + -0.041014*tf.T953 + -1.5*tf.lnT9)
+    # ac12n
+    rate += np.exp(  -8.22898 + -13.3317*tf.T913i + 13.2237*tf.T913
+                  + -9.06339*tf.T9 + -0.666667*tf.lnT9)
+
+    rate_eval.n_He4_He4__Be9 = rate
+
 @numba.njit()
 def He4_He4_He4__C12(rate_eval, tf):
     # 3 He4 --> C12
@@ -7739,6 +8188,23 @@ def p_He4_He4__d_Be7(rate_eval, tf):
 
     rate_eval.p_He4_He4__d_Be7 = rate
 
+@numba.njit()
+def d_He4_He4__p_Be9(rate_eval, tf):
+    # d + He4 + He4 --> p + Be9
+    rate = 0.0
+
+    # cf88n
+    rate += np.exp(  2.83369 + -7.55453*tf.T9i + -10.359*tf.T913i + 0.103955*tf.T913
+                  + 4.4262*tf.T9 + -5.95664*tf.T953 + -2.16667*tf.lnT9)
+    # cf88r
+    rate += np.exp(  -2.68071 + -13.3545*tf.T9i
+                  + -2.25*tf.lnT9)
+    # cf88r
+    rate += np.exp(  -3.06461 + -10.6005*tf.T9i
+                  + -2.5*tf.lnT9)
+
+    rate_eval.d_He4_He4__p_Be9 = rate
+
 @numba.njit()
 def He4_He4_He4__n_C11(rate_eval, tf):
     # 3 He4 --> n + C11
@@ -7789,6 +8255,20 @@ def n_p_He4_He4__t_Be7(rate_eval, tf):
 
     rate_eval.n_p_He4_He4__t_Be7 = rate
 
+@numba.njit()
+def n_p_He4_He4__p_Be9(rate_eval, tf):
+    # n + p + He4 + He4 --> p + Be9
+    rate = 0.0
+
+    # cf88r
+    rate += np.exp(  -25.452 + -8.47112*tf.T9i + 1.40505e-06*tf.T913i + -1.47128e-06*tf.T913
+                  + 6.89313e-08*tf.T9 + -3.55179e-09*tf.T953 + -4.5*tf.lnT9)
+    # cf88n
+    rate += np.exp(  -31.5916 + -3.22032*tf.T9i + -0.634849*tf.T913i + 4.82033*tf.T913
+                  + -0.257317*tf.T9 + 0.0134206*tf.T953 + -4.08885*tf.lnT9)
+
+    rate_eval.n_p_He4_He4__p_Be9 = rate
+
 @numba.njit()
 def p_p_He4_He4__He3_Be7(rate_eval, tf):
     # p + p + He4 + He4 --> He3 + Be7
@@ -7801,7 +8281,6 @@ def p_p_He4_He4__He3_Be7(rate_eval, tf):
     rate_eval.p_p_He4_He4__He3_Be7 = rate
 
 def rhs(t, Y, rho, T, screen_func=None):
-    print(f"TIME - {t:0.3E}")
     return rhs_eq(t, Y, rho, T, screen_func)
 
 @numba.njit()
@@ -7814,7 +8293,9 @@ def rhs_eq(t, Y, rho, T, screen_func):
     n__p__weak__wc12(rate_eval, tf)
     t__He3__weak__wc12(rate_eval, tf)
     He3__t__weak__electron_capture(rate_eval, tf)
+    He6__Li6__weak__wc12(rate_eval, tf)
     Be7__Li7__weak__electron_capture(rate_eval, tf)
+    C14__N14__weak__wc12(rate_eval, tf)
     N13__C13__weak__wc12(rate_eval, tf)
     O14__N14__weak__wc12(rate_eval, tf)
     O15__N15__weak__wc12(rate_eval, tf)
@@ -7853,10 +8334,12 @@ def rhs_eq(t, Y, rho, T, screen_func):
     C11__He4_Be7(rate_eval, tf)
     C12__n_C11(rate_eval, tf)
     C13__n_C12(rate_eval, tf)
+    C14__n_C13(rate_eval, tf)
     N13__p_C12(rate_eval, tf)
     N14__n_N13(rate_eval, tf)
     N14__p_C13(rate_eval, tf)
     N15__n_N14(rate_eval, tf)
+    N15__p_C14(rate_eval, tf)
     O14__p_N13(rate_eval, tf)
     O15__n_O14(rate_eval, tf)
     O15__p_N14(rate_eval, tf)
@@ -7865,6 +8348,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
     O16__He4_C12(rate_eval, tf)
     O17__n_O16(rate_eval, tf)
     O18__n_O17(rate_eval, tf)
+    O18__He4_C14(rate_eval, tf)
     F15__p_O14(rate_eval, tf)
     F16__n_F15(rate_eval, tf)
     F16__p_O15(rate_eval, tf)
@@ -7949,7 +8433,9 @@ def rhs_eq(t, Y, rho, T, screen_func):
     Cl35__He4_P31(rate_eval, tf)
     Ar35__p_Cl34(rate_eval, tf)
     Ar35__He4_S31(rate_eval, tf)
+    He6__n_n_He4(rate_eval, tf)
     Li6__n_p_He4(rate_eval, tf)
+    Be9__n_He4_He4(rate_eval, tf)
     C12__He4_He4_He4(rate_eval, tf)
     n_p__d(rate_eval, tf)
     p_p__d__weak__bet_pos_(rate_eval, tf)
@@ -7971,7 +8457,10 @@ def rhs_eq(t, Y, rho, T, screen_func):
     n_C12__C13(rate_eval, tf)
     p_C12__N13(rate_eval, tf)
     He4_C12__O16(rate_eval, tf)
+    n_C13__C14(rate_eval, tf)
     p_C13__N14(rate_eval, tf)
+    p_C14__N15(rate_eval, tf)
+    He4_C14__O18(rate_eval, tf)
     n_N13__N14(rate_eval, tf)
     p_N13__O14(rate_eval, tf)
     n_N14__N15(rate_eval, tf)
@@ -8083,16 +8572,22 @@ def rhs_eq(t, Y, rho, T, screen_func):
     p_Li6__He4_He3(rate_eval, tf)
     d_Li6__n_Be7(rate_eval, tf)
     d_Li6__p_Li7(rate_eval, tf)
+    He4_Li6__p_Be9(rate_eval, tf)
     p_Li7__n_Be7(rate_eval, tf)
     p_Li7__d_Li6(rate_eval, tf)
     p_Li7__He4_He4(rate_eval, tf)
+    t_Li7__n_Be9(rate_eval, tf)
     n_Be7__p_Li7(rate_eval, tf)
     n_Be7__d_Li6(rate_eval, tf)
     n_Be7__He4_He4(rate_eval, tf)
+    n_Be9__t_Li7(rate_eval, tf)
+    p_Be9__He4_Li6(rate_eval, tf)
+    He4_Be9__n_C12(rate_eval, tf)
     He4_B8__p_C11(rate_eval, tf)
     p_C11__He4_B8(rate_eval, tf)
     He4_C11__n_O14(rate_eval, tf)
     He4_C11__p_N14(rate_eval, tf)
+    n_C12__He4_Be9(rate_eval, tf)
     He4_C12__n_O15(rate_eval, tf)
     He4_C12__p_N15(rate_eval, tf)
     C12_C12__n_Mg23(rate_eval, tf)
@@ -8101,14 +8596,19 @@ def rhs_eq(t, Y, rho, T, screen_func):
     p_C13__n_N13(rate_eval, tf)
     d_C13__n_N14(rate_eval, tf)
     He4_C13__n_O16(rate_eval, tf)
+    p_C14__n_N14(rate_eval, tf)
+    d_C14__n_N15(rate_eval, tf)
+    He4_C14__n_O17(rate_eval, tf)
     n_N13__p_C13(rate_eval, tf)
     He4_N13__n_F16(rate_eval, tf)
     He4_N13__p_O16(rate_eval, tf)
+    n_N14__p_C14(rate_eval, tf)
     n_N14__d_C13(rate_eval, tf)
     p_N14__n_O14(rate_eval, tf)
     p_N14__He4_C11(rate_eval, tf)
     He4_N14__n_F17(rate_eval, tf)
     He4_N14__p_O17(rate_eval, tf)
+    n_N15__d_C14(rate_eval, tf)
     p_N15__n_O15(rate_eval, tf)
     p_N15__He4_C12(rate_eval, tf)
     He4_N15__n_F18(rate_eval, tf)
@@ -8131,6 +8631,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
     O16_O16__n_S31(rate_eval, tf)
     O16_O16__p_P31(rate_eval, tf)
     O16_O16__He4_Si28(rate_eval, tf)
+    n_O17__He4_C14(rate_eval, tf)
     p_O17__n_F17(rate_eval, tf)
     p_O17__He4_N14(rate_eval, tf)
     He4_O17__n_Ne20(rate_eval, tf)
@@ -8297,13 +8798,17 @@ def rhs_eq(t, Y, rho, T, screen_func):
     He3_He3__p_p_He4(rate_eval, tf)
     d_Li7__n_He4_He4(rate_eval, tf)
     d_Be7__p_He4_He4(rate_eval, tf)
+    p_Be9__d_He4_He4(rate_eval, tf)
     n_B8__p_He4_He4(rate_eval, tf)
     n_C11__He4_He4_He4(rate_eval, tf)
     t_Li7__n_n_He4_He4(rate_eval, tf)
     He3_Li7__n_p_He4_He4(rate_eval, tf)
     t_Be7__n_p_He4_He4(rate_eval, tf)
     He3_Be7__p_p_He4_He4(rate_eval, tf)
+    p_Be9__n_p_He4_He4(rate_eval, tf)
+    n_n_He4__He6(rate_eval, tf)
     n_p_He4__Li6(rate_eval, tf)
+    n_He4_He4__Be9(rate_eval, tf)
     He4_He4_He4__C12(rate_eval, tf)
     n_p_p__p_d(rate_eval, tf)
     n_n_He4__t_t(rate_eval, tf)
@@ -8312,10 +8817,12 @@ def rhs_eq(t, Y, rho, T, screen_func):
     n_He4_He4__d_Li7(rate_eval, tf)
     p_He4_He4__n_B8(rate_eval, tf)
     p_He4_He4__d_Be7(rate_eval, tf)
+    d_He4_He4__p_Be9(rate_eval, tf)
     He4_He4_He4__n_C11(rate_eval, tf)
     n_n_He4_He4__t_Li7(rate_eval, tf)
     n_p_He4_He4__He3_Li7(rate_eval, tf)
     n_p_He4_He4__t_Be7(rate_eval, tf)
+    n_p_He4_He4__p_Be9(rate_eval, tf)
     p_p_He4_He4__He3_Be7(rate_eval, tf)
 
     if screen_func is not None:
@@ -8392,6 +8899,16 @@ def rhs_eq(t, Y, rho, T, screen_func):
         rate_eval.p_C13__N14 *= scor
         rate_eval.p_C13__n_N13 *= scor
 
+        scn_fac = ScreenFactors(1, 1, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.p_C14__N15 *= scor
+        rate_eval.p_C14__n_N14 *= scor
+
+        scn_fac = ScreenFactors(2, 4, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_C14__O18 *= scor
+        rate_eval.He4_C14__n_O17 *= scor
+
         scn_fac = ScreenFactors(1, 1, 7, 13)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.p_N13__O14 *= scor
@@ -8755,6 +9272,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_He4__n_Be7 *= scor
         rate_eval.He4_He4__p_Li7 *= scor
+        rate_eval.n_He4_He4__Be9 *= scor
         rate_eval.n_He4_He4__d_Li7 *= scor
         rate_eval.n_n_He4_He4__t_Li7 *= scor
 
@@ -8763,12 +9281,31 @@ def rhs_eq(t, Y, rho, T, screen_func):
         rate_eval.d_Li6__n_Be7 *= scor
         rate_eval.d_Li6__p_Li7 *= scor
 
+        scn_fac = ScreenFactors(2, 4, 3, 6)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_Li6__p_Be9 *= scor
+
         scn_fac = ScreenFactors(1, 1, 3, 7)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.p_Li7__n_Be7 *= scor
         rate_eval.p_Li7__d_Li6 *= scor
         rate_eval.p_Li7__He4_He4 *= scor
 
+        scn_fac = ScreenFactors(1, 3, 3, 7)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.t_Li7__n_Be9 *= scor
+        rate_eval.t_Li7__n_n_He4_He4 *= scor
+
+        scn_fac = ScreenFactors(1, 1, 4, 9)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.p_Be9__He4_Li6 *= scor
+        rate_eval.p_Be9__d_He4_He4 *= scor
+        rate_eval.p_Be9__n_p_He4_He4 *= scor
+
+        scn_fac = ScreenFactors(2, 4, 4, 9)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_Be9__n_C12 *= scor
+
         scn_fac = ScreenFactors(2, 4, 5, 8)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_B8__p_C11 *= scor
@@ -8796,6 +9333,10 @@ def rhs_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_C13__n_O16 *= scor
 
+        scn_fac = ScreenFactors(1, 2, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.d_C14__n_N15 *= scor
+
         scn_fac = ScreenFactors(2, 4, 7, 13)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_N13__n_F16 *= scor
@@ -8914,10 +9455,6 @@ def rhs_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.d_Be7__p_He4_He4 *= scor
 
-        scn_fac = ScreenFactors(1, 3, 3, 7)
-        scor = screen_func(plasma_state, scn_fac)
-        rate_eval.t_Li7__n_n_He4_He4 *= scor
-
         scn_fac = ScreenFactors(2, 3, 3, 7)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He3_Li7__n_p_He4_He4 *= scor
@@ -8947,6 +9484,11 @@ def rhs_eq(t, Y, rho, T, screen_func):
         rate_eval.p_He4_He4__d_Be7 *= scor
         rate_eval.n_p_He4_He4__He3_Li7 *= scor
         rate_eval.n_p_He4_He4__t_Be7 *= scor
+        rate_eval.n_p_He4_He4__p_Be9 *= scor
+
+        scn_fac = ScreenFactors(1, 2, 2, 4)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.d_He4_He4__p_Be9 *= scor
 
     dYdt = np.zeros((nnuc), dtype=np.float64)
 
@@ -8958,6 +9500,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jn]*Y[jli6]*rate_eval.n_Li6__Li7 +Y[jli7]*rate_eval.Li7__n_Li6 ) +
           ( -rho*Y[jn]*Y[jc11]*rate_eval.n_C11__C12 +Y[jc12]*rate_eval.C12__n_C11 ) +
           ( -rho*Y[jn]*Y[jc12]*rate_eval.n_C12__C13 +Y[jc13]*rate_eval.C13__n_C12 ) +
+          ( -rho*Y[jn]*Y[jc13]*rate_eval.n_C13__C14 +Y[jc14]*rate_eval.C14__n_C13 ) +
           ( -rho*Y[jn]*Y[jn13]*rate_eval.n_N13__N14 +Y[jn14]*rate_eval.N14__n_N13 ) +
           ( -rho*Y[jn]*Y[jn14]*rate_eval.n_N14__N15 +Y[jn15]*rate_eval.N15__n_N14 ) +
           ( -rho*Y[jn]*Y[jo14]*rate_eval.n_O14__O15 +Y[jo15]*rate_eval.O15__n_O14 ) +
@@ -8993,16 +9536,21 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jn]*Y[jhe3]*rate_eval.n_He3__p_t +rho*Y[jp]*Y[jt]*rate_eval.p_t__n_He3 ) +
           ( -rho*Y[jn]*Y[jli6]*rate_eval.n_Li6__He4_t +rho*Y[jt]*Y[jhe4]*rate_eval.He4_t__n_Li6 ) +
           ( +rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__n_Be7 -rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__d_Li6 ) +
+          ( +rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_Be9 -rho*Y[jn]*Y[jbe9]*rate_eval.n_Be9__t_Li7 ) +
           ( -rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__p_Li7 +rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__n_Be7 ) +
           ( -rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__He4_He4 +5.00000000000000e-01*rho*Y[jhe4]**2*rate_eval.He4_He4__n_Be7 ) +
+          ( +rho*Y[jhe4]*Y[jbe9]*rate_eval.He4_Be9__n_C12 -rho*Y[jn]*Y[jc12]*rate_eval.n_C12__He4_Be9 ) +
           ( +rho*Y[jd]*Y[jc13]*rate_eval.d_C13__n_N14 -rho*Y[jn]*Y[jn14]*rate_eval.n_N14__d_C13 ) +
           ( +rho*Y[jhe4]*Y[jc13]*rate_eval.He4_C13__n_O16 -rho*Y[jn]*Y[jo16]*rate_eval.n_O16__He4_C13 ) +
+          ( +rho*Y[jd]*Y[jc14]*rate_eval.d_C14__n_N15 -rho*Y[jn]*Y[jn15]*rate_eval.n_N15__d_C14 ) +
           ( -rho*Y[jn]*Y[jn13]*rate_eval.n_N13__p_C13 +rho*Y[jp]*Y[jc13]*rate_eval.p_C13__n_N13 ) +
+          ( -rho*Y[jn]*Y[jn14]*rate_eval.n_N14__p_C14 +rho*Y[jp]*Y[jc14]*rate_eval.p_C14__n_N14 ) +
           ( -rho*Y[jn]*Y[jo14]*rate_eval.n_O14__p_N14 +rho*Y[jp]*Y[jn14]*rate_eval.p_N14__n_O14 ) +
           ( -rho*Y[jn]*Y[jo14]*rate_eval.n_O14__He4_C11 +rho*Y[jhe4]*Y[jc11]*rate_eval.He4_C11__n_O14 ) +
           ( -rho*Y[jn]*Y[jo15]*rate_eval.n_O15__p_N15 +rho*Y[jp]*Y[jn15]*rate_eval.p_N15__n_O15 ) +
           ( -rho*Y[jn]*Y[jo15]*rate_eval.n_O15__He4_C12 +rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__n_O15 ) +
           ( +5.00000000000000e-01*rho*Y[jo16]**2*rate_eval.O16_O16__n_S31 -rho*Y[jn]*Y[js31]*rate_eval.n_S31__O16_O16 ) +
+          ( -rho*Y[jn]*Y[jo17]*rate_eval.n_O17__He4_C14 +rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__n_O17 ) +
           ( +rho*Y[jhe4]*Y[jo17]*rate_eval.He4_O17__n_Ne20 -rho*Y[jn]*Y[jne20]*rate_eval.n_Ne20__He4_O17 ) +
           ( -rho*Y[jn]*Y[jf15]*rate_eval.n_F15__p_O15 +rho*Y[jp]*Y[jo15]*rate_eval.p_O15__n_F15 ) +
           ( -rho*Y[jn]*Y[jf16]*rate_eval.n_F16__p_O16 +rho*Y[jp]*Y[jo16]*rate_eval.p_O16__n_F16 ) +
@@ -9061,8 +9609,11 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( + 2*rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4 + -2*2.50000000000000e-01*rho**3*Y[jn]**2*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7 ) +
           ( +rho*Y[jhe3]*Y[jli7]*rate_eval.He3_Li7__n_p_He4_He4 -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7 ) +
           ( +rho*Y[jt]*Y[jbe7]*rate_eval.t_Be7__n_p_He4_He4 -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7 ) +
+          ( + -2*5.00000000000000e-01*rho**2*Y[jn]**2*Y[jhe4]*rate_eval.n_n_He4__He6 + 2*Y[jhe6]*rate_eval.He6__n_n_He4 ) +
           ( -rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6 +Y[jli6]*rate_eval.Li6__n_p_He4 ) +
-          ( -5.00000000000000e-01*rho**2*Y[jn]*Y[jp]**2*rate_eval.n_p_p__p_d +rho*Y[jp]*Y[jd]*rate_eval.p_d__n_p_p )
+          ( -5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__Be9 +Y[jbe9]*rate_eval.Be9__n_He4_He4 ) +
+          ( -5.00000000000000e-01*rho**2*Y[jn]*Y[jp]**2*rate_eval.n_p_p__p_d +rho*Y[jp]*Y[jd]*rate_eval.p_d__n_p_p ) +
+          ( -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9 +rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4 )
        )
 
     dYdt[jp] = (
@@ -9080,6 +9631,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jp]*Y[jbe7]*rate_eval.p_Be7__B8 +Y[jb8]*rate_eval.B8__p_Be7 ) +
           ( -rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13 +Y[jn13]*rate_eval.N13__p_C12 ) +
           ( -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14 +Y[jn14]*rate_eval.N14__p_C13 ) +
+          ( -rho*Y[jp]*Y[jc14]*rate_eval.p_C14__N15 +Y[jn15]*rate_eval.N15__p_C14 ) +
           ( -rho*Y[jp]*Y[jn13]*rate_eval.p_N13__O14 +Y[jo14]*rate_eval.O14__p_N13 ) +
           ( -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__O15 +Y[jo15]*rate_eval.O15__p_N14 ) +
           ( -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__O16 +Y[jo16]*rate_eval.O16__p_N15 ) +
@@ -9116,11 +9668,13 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__p_Li7 -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__d_Li6 ) +
           ( -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__He4_He4 +5.00000000000000e-01*rho*Y[jhe4]**2*rate_eval.He4_He4__p_Li7 ) +
           ( +rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__p_Li7 -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__n_Be7 ) +
+          ( -rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__He4_Li6 +rho*Y[jhe4]*Y[jli6]*rate_eval.He4_Li6__p_Be9 ) +
           ( +rho*Y[jhe4]*Y[jb8]*rate_eval.He4_B8__p_C11 -rho*Y[jp]*Y[jc11]*rate_eval.p_C11__He4_B8 ) +
           ( +rho*Y[jhe4]*Y[jc11]*rate_eval.He4_C11__p_N14 -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__He4_C11 ) +
           ( +5.00000000000000e-01*rho*Y[jc12]**2*rate_eval.C12_C12__p_Na23 -rho*Y[jp]*Y[jna23]*rate_eval.p_Na23__C12_C12 ) +
           ( +rho*Y[jn]*Y[jn13]*rate_eval.n_N13__p_C13 -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__n_N13 ) +
           ( +rho*Y[jhe4]*Y[jn13]*rate_eval.He4_N13__p_O16 -rho*Y[jp]*Y[jo16]*rate_eval.p_O16__He4_N13 ) +
+          ( +rho*Y[jn]*Y[jn14]*rate_eval.n_N14__p_C14 -rho*Y[jp]*Y[jc14]*rate_eval.p_C14__n_N14 ) +
           ( -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12 +rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__p_N15 ) +
           ( +rho*Y[jn]*Y[jo14]*rate_eval.n_O14__p_N14 -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__n_O14 ) +
           ( +rho*Y[jhe4]*Y[jo14]*rate_eval.He4_O14__p_F17 -rho*Y[jp]*Y[jf17]*rate_eval.p_F17__He4_O14 ) +
@@ -9182,12 +9736,14 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jt]*Y[jhe3]*rate_eval.t_He3__n_p_He4 -rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3 ) +
           ( + 2*5.00000000000000e-01*rho*Y[jhe3]**2*rate_eval.He3_He3__p_p_He4 + -2*5.00000000000000e-01*rho**2*Y[jp]**2*Y[jhe4]*rate_eval.p_p_He4__He3_He3 ) +
           ( +rho*Y[jd]*Y[jbe7]*rate_eval.d_Be7__p_He4_He4 -5.00000000000000e-01*rho**2*Y[jp]*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7 ) +
+          ( -rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__d_He4_He4 +5.00000000000000e-01*rho**2*Y[jd]*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9 ) +
           ( +rho*Y[jn]*Y[jb8]*rate_eval.n_B8__p_He4_He4 -5.00000000000000e-01*rho**2*Y[jp]*Y[jhe4]**2*rate_eval.p_He4_He4__n_B8 ) +
           ( +rho*Y[jhe3]*Y[jli7]*rate_eval.He3_Li7__n_p_He4_He4 -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7 ) +
           ( +rho*Y[jt]*Y[jbe7]*rate_eval.t_Be7__n_p_He4_He4 -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7 ) +
           ( + 2*rho*Y[jhe3]*Y[jbe7]*rate_eval.He3_Be7__p_p_He4_He4 + -2*2.50000000000000e-01*rho**3*Y[jp]**2*Y[jhe4]**2*rate_eval.p_p_He4_He4__He3_Be7 ) +
           ( -rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6 +Y[jli6]*rate_eval.Li6__n_p_He4 ) +
-          ( -5.00000000000000e-01*rho**2*Y[jn]*Y[jp]**2*rate_eval.n_p_p__p_d +rho*Y[jp]*Y[jd]*rate_eval.p_d__n_p_p )
+          ( -5.00000000000000e-01*rho**2*Y[jn]*Y[jp]**2*rate_eval.n_p_p__p_d +rho*Y[jp]*Y[jd]*rate_eval.p_d__n_p_p ) +
+          ( + 0*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9 + 0*rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4 )
        )
 
     dYdt[jd] = (
@@ -9206,8 +9762,10 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__n_Be7 +rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__d_Li6 ) +
           ( -rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__p_Li7 +rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__d_Li6 ) +
           ( -rho*Y[jd]*Y[jc13]*rate_eval.d_C13__n_N14 +rho*Y[jn]*Y[jn14]*rate_eval.n_N14__d_C13 ) +
+          ( -rho*Y[jd]*Y[jc14]*rate_eval.d_C14__n_N15 +rho*Y[jn]*Y[jn15]*rate_eval.n_N15__d_C14 ) +
           ( -rho*Y[jd]*Y[jli7]*rate_eval.d_Li7__n_He4_He4 +5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7 ) +
           ( -rho*Y[jd]*Y[jbe7]*rate_eval.d_Be7__p_He4_He4 +5.00000000000000e-01*rho**2*Y[jp]*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7 ) +
+          ( +rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__d_He4_He4 -5.00000000000000e-01*rho**2*Y[jd]*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9 ) +
           ( +5.00000000000000e-01*rho**2*Y[jn]*Y[jp]**2*rate_eval.n_p_p__p_d -rho*Y[jp]*Y[jd]*rate_eval.p_d__n_p_p )
        )
 
@@ -9221,6 +9779,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jn]*Y[jhe3]*rate_eval.n_He3__p_t -rho*Y[jp]*Y[jt]*rate_eval.p_t__n_He3 ) +
           ( -rho*Y[jt]*Y[jhe3]*rate_eval.t_He3__d_He4 +rho*Y[jd]*Y[jhe4]*rate_eval.d_He4__t_He3 ) +
           ( +rho*Y[jn]*Y[jli6]*rate_eval.n_Li6__He4_t -rho*Y[jt]*Y[jhe4]*rate_eval.He4_t__n_Li6 ) +
+          ( -rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_Be9 +rho*Y[jn]*Y[jbe9]*rate_eval.n_Be9__t_Li7 ) +
           ( + -2*5.00000000000000e-01*rho*Y[jt]**2*rate_eval.t_t__n_n_He4 + 2*5.00000000000000e-01*rho**2*Y[jn]**2*Y[jhe4]*rate_eval.n_n_He4__t_t ) +
           ( -rho*Y[jt]*Y[jhe3]*rate_eval.t_He3__n_p_He4 +rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3 ) +
           ( -rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4 +2.50000000000000e-01*rho**3*Y[jn]**2*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7 ) +
@@ -9256,6 +9815,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jhe3]*Y[jhe4]*rate_eval.He4_He3__Be7 +Y[jbe7]*rate_eval.Be7__He4_He3 ) +
           ( -rho*Y[jhe4]*Y[jbe7]*rate_eval.He4_Be7__C11 +Y[jc11]*rate_eval.C11__He4_Be7 ) +
           ( -rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__O16 +Y[jo16]*rate_eval.O16__He4_C12 ) +
+          ( -rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__O18 +Y[jo18]*rate_eval.O18__He4_C14 ) +
           ( -rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__F18 +Y[jf18]*rate_eval.F18__He4_N14 ) +
           ( -rho*Y[jhe4]*Y[jn15]*rate_eval.He4_N15__F19 +Y[jf19]*rate_eval.F19__He4_N15 ) +
           ( -rho*Y[jhe4]*Y[jo14]*rate_eval.He4_O14__Ne18 +Y[jne18]*rate_eval.Ne18__He4_O14 ) +
@@ -9293,6 +9853,8 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jp]*Y[jli6]*rate_eval.p_Li6__He4_He3 -rho*Y[jhe3]*Y[jhe4]*rate_eval.He4_He3__p_Li6 ) +
           ( + 2*rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__He4_He4 + -2*5.00000000000000e-01*rho*Y[jhe4]**2*rate_eval.He4_He4__p_Li7 ) +
           ( + 2*rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__He4_He4 + -2*5.00000000000000e-01*rho*Y[jhe4]**2*rate_eval.He4_He4__n_Be7 ) +
+          ( +rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__He4_Li6 -rho*Y[jhe4]*Y[jli6]*rate_eval.He4_Li6__p_Be9 ) +
+          ( -rho*Y[jhe4]*Y[jbe9]*rate_eval.He4_Be9__n_C12 +rho*Y[jn]*Y[jc12]*rate_eval.n_C12__He4_Be9 ) +
           ( -rho*Y[jhe4]*Y[jb8]*rate_eval.He4_B8__p_C11 +rho*Y[jp]*Y[jc11]*rate_eval.p_C11__He4_B8 ) +
           ( -rho*Y[jhe4]*Y[jc11]*rate_eval.He4_C11__p_N14 +rho*Y[jp]*Y[jn14]*rate_eval.p_N14__He4_C11 ) +
           ( +5.00000000000000e-01*rho*Y[jc12]**2*rate_eval.C12_C12__He4_Ne20 -rho*Y[jhe4]*Y[jne20]*rate_eval.He4_Ne20__C12_C12 ) +
@@ -9304,6 +9866,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jn]*Y[jo15]*rate_eval.n_O15__He4_C12 -rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__n_O15 ) +
           ( +rho*Y[jc12]*Y[jo16]*rate_eval.C12_O16__He4_Mg24 -rho*Y[jhe4]*Y[jmg24]*rate_eval.He4_Mg24__C12_O16 ) +
           ( +5.00000000000000e-01*rho*Y[jo16]**2*rate_eval.O16_O16__He4_Si28 -rho*Y[jhe4]*Y[jsi28]*rate_eval.He4_Si28__O16_O16 ) +
+          ( +rho*Y[jn]*Y[jo17]*rate_eval.n_O17__He4_C14 -rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__n_O17 ) +
           ( +rho*Y[jp]*Y[jo17]*rate_eval.p_O17__He4_N14 -rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__p_O17 ) +
           ( -rho*Y[jhe4]*Y[jo17]*rate_eval.He4_O17__n_Ne20 +rho*Y[jn]*Y[jne20]*rate_eval.n_Ne20__He4_O17 ) +
           ( +rho*Y[jp]*Y[jo18]*rate_eval.p_O18__He4_N15 -rho*Y[jhe4]*Y[jn15]*rate_eval.He4_N15__p_O18 ) +
@@ -9367,17 +9930,27 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +5.00000000000000e-01*rho*Y[jhe3]**2*rate_eval.He3_He3__p_p_He4 -5.00000000000000e-01*rho**2*Y[jp]**2*Y[jhe4]*rate_eval.p_p_He4__He3_He3 ) +
           ( + 2*rho*Y[jd]*Y[jli7]*rate_eval.d_Li7__n_He4_He4 + -2*5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7 ) +
           ( + 2*rho*Y[jd]*Y[jbe7]*rate_eval.d_Be7__p_He4_He4 + -2*5.00000000000000e-01*rho**2*Y[jp]*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7 ) +
+          ( + 2*rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__d_He4_He4 + -2*5.00000000000000e-01*rho**2*Y[jd]*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9 ) +
           ( + 2*rho*Y[jn]*Y[jb8]*rate_eval.n_B8__p_He4_He4 + -2*5.00000000000000e-01*rho**2*Y[jp]*Y[jhe4]**2*rate_eval.p_He4_He4__n_B8 ) +
           ( + 3*rho*Y[jn]*Y[jc11]*rate_eval.n_C11__He4_He4_He4 + -3*1.66666666666667e-01*rho**2*Y[jhe4]**3*rate_eval.He4_He4_He4__n_C11 ) +
           ( + 2*rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4 + -2*2.50000000000000e-01*rho**3*Y[jn]**2*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7 ) +
           ( + 2*rho*Y[jhe3]*Y[jli7]*rate_eval.He3_Li7__n_p_He4_He4 + -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7 ) +
           ( + 2*rho*Y[jt]*Y[jbe7]*rate_eval.t_Be7__n_p_He4_He4 + -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7 ) +
           ( + 2*rho*Y[jhe3]*Y[jbe7]*rate_eval.He3_Be7__p_p_He4_He4 + -2*2.50000000000000e-01*rho**3*Y[jp]**2*Y[jhe4]**2*rate_eval.p_p_He4_He4__He3_Be7 ) +
+          ( -5.00000000000000e-01*rho**2*Y[jn]**2*Y[jhe4]*rate_eval.n_n_He4__He6 +Y[jhe6]*rate_eval.He6__n_n_He4 ) +
           ( -rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6 +Y[jli6]*rate_eval.Li6__n_p_He4 ) +
-          ( + -3*1.66666666666667e-01*rho**2*Y[jhe4]**3*rate_eval.He4_He4_He4__C12 + 3*Y[jc12]*rate_eval.C12__He4_He4_He4 )
+          ( + -2*5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__Be9 + 2*Y[jbe9]*rate_eval.Be9__n_He4_He4 ) +
+          ( + -3*1.66666666666667e-01*rho**2*Y[jhe4]**3*rate_eval.He4_He4_He4__C12 + 3*Y[jc12]*rate_eval.C12__He4_He4_He4 ) +
+          ( + -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9 + 2*rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4 )
+       )
+
+    dYdt[jhe6] = (
+          -Y[jhe6]*rate_eval.He6__Li6__weak__wc12  +
+          ( +5.00000000000000e-01*rho**2*Y[jn]**2*Y[jhe4]*rate_eval.n_n_He4__He6 -Y[jhe6]*rate_eval.He6__n_n_He4 )
        )
 
     dYdt[jli6] = (
+          +Y[jhe6]*rate_eval.He6__Li6__weak__wc12  +
           ( +rho*Y[jd]*Y[jhe4]*rate_eval.He4_d__Li6 -Y[jli6]*rate_eval.Li6__He4_d ) +
           ( -rho*Y[jn]*Y[jli6]*rate_eval.n_Li6__Li7 +Y[jli7]*rate_eval.Li7__n_Li6 ) +
           ( -rho*Y[jp]*Y[jli6]*rate_eval.p_Li6__Be7 +Y[jbe7]*rate_eval.Be7__p_Li6 ) +
@@ -9385,6 +9958,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jp]*Y[jli6]*rate_eval.p_Li6__He4_He3 +rho*Y[jhe3]*Y[jhe4]*rate_eval.He4_He3__p_Li6 ) +
           ( -rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__n_Be7 +rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__d_Li6 ) +
           ( -rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__p_Li7 +rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__d_Li6 ) +
+          ( +rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__He4_Li6 -rho*Y[jhe4]*Y[jli6]*rate_eval.He4_Li6__p_Be9 ) +
           ( +rho**2*Y[jn]*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6 -Y[jli6]*rate_eval.Li6__n_p_He4 )
        )
 
@@ -9394,6 +9968,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( +rho*Y[jn]*Y[jli6]*rate_eval.n_Li6__Li7 -Y[jli7]*rate_eval.Li7__n_Li6 ) +
           ( +rho*Y[jd]*Y[jli6]*rate_eval.d_Li6__p_Li7 -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__d_Li6 ) +
           ( -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__He4_He4 +5.00000000000000e-01*rho*Y[jhe4]**2*rate_eval.He4_He4__p_Li7 ) +
+          ( -rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_Be9 +rho*Y[jn]*Y[jbe9]*rate_eval.n_Be9__t_Li7 ) +
           ( +rho*Y[jn]*Y[jbe7]*rate_eval.n_Be7__p_Li7 -rho*Y[jp]*Y[jli7]*rate_eval.p_Li7__n_Be7 ) +
           ( -rho*Y[jd]*Y[jli7]*rate_eval.d_Li7__n_He4_He4 +5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7 ) +
           ( -rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4 +2.50000000000000e-01*rho**3*Y[jn]**2*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7 ) +
@@ -9414,6 +9989,15 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jhe3]*Y[jbe7]*rate_eval.He3_Be7__p_p_He4_He4 +2.50000000000000e-01*rho**3*Y[jp]**2*Y[jhe4]**2*rate_eval.p_p_He4_He4__He3_Be7 )
        )
 
+    dYdt[jbe9] = (
+          ( +rho*Y[jt]*Y[jli7]*rate_eval.t_Li7__n_Be9 -rho*Y[jn]*Y[jbe9]*rate_eval.n_Be9__t_Li7 ) +
+          ( -rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__He4_Li6 +rho*Y[jhe4]*Y[jli6]*rate_eval.He4_Li6__p_Be9 ) +
+          ( -rho*Y[jhe4]*Y[jbe9]*rate_eval.He4_Be9__n_C12 +rho*Y[jn]*Y[jc12]*rate_eval.n_C12__He4_Be9 ) +
+          ( -rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__d_He4_He4 +5.00000000000000e-01*rho**2*Y[jd]*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9 ) +
+          ( +5.00000000000000e-01*rho**2*Y[jn]*Y[jhe4]**2*rate_eval.n_He4_He4__Be9 -Y[jbe9]*rate_eval.Be9__n_He4_He4 ) +
+          ( +5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9 -rho*Y[jp]*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4 )
+       )
+
     dYdt[jb8] = (
           -Y[jb8]*rate_eval.B8__He4_He4__weak__wc12  +
           ( +rho*Y[jp]*Y[jbe7]*rate_eval.p_Be7__B8 -Y[jb8]*rate_eval.B8__p_Be7 ) +
@@ -9435,6 +10019,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jn]*Y[jc12]*rate_eval.n_C12__C13 +Y[jc13]*rate_eval.C13__n_C12 ) +
           ( -rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13 +Y[jn13]*rate_eval.N13__p_C12 ) +
           ( -rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__O16 +Y[jo16]*rate_eval.O16__He4_C12 ) +
+          ( +rho*Y[jhe4]*Y[jbe9]*rate_eval.He4_Be9__n_C12 -rho*Y[jn]*Y[jc12]*rate_eval.n_C12__He4_Be9 ) +
           ( + -2*5.00000000000000e-01*rho*Y[jc12]**2*rate_eval.C12_C12__p_Na23 + 2*rho*Y[jp]*Y[jna23]*rate_eval.p_Na23__C12_C12 ) +
           ( + -2*5.00000000000000e-01*rho*Y[jc12]**2*rate_eval.C12_C12__He4_Ne20 + 2*rho*Y[jhe4]*Y[jne20]*rate_eval.He4_Ne20__C12_C12 ) +
           ( +rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12 -rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__p_N15 ) +
@@ -9451,12 +10036,23 @@ def rhs_eq(t, Y, rho, T, screen_func):
     dYdt[jc13] = (
           +Y[jn13]*rate_eval.N13__C13__weak__wc12  +
           ( +rho*Y[jn]*Y[jc12]*rate_eval.n_C12__C13 -Y[jc13]*rate_eval.C13__n_C12 ) +
+          ( -rho*Y[jn]*Y[jc13]*rate_eval.n_C13__C14 +Y[jc14]*rate_eval.C14__n_C13 ) +
           ( -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14 +Y[jn14]*rate_eval.N14__p_C13 ) +
           ( -rho*Y[jd]*Y[jc13]*rate_eval.d_C13__n_N14 +rho*Y[jn]*Y[jn14]*rate_eval.n_N14__d_C13 ) +
           ( -rho*Y[jhe4]*Y[jc13]*rate_eval.He4_C13__n_O16 +rho*Y[jn]*Y[jo16]*rate_eval.n_O16__He4_C13 ) +
           ( +rho*Y[jn]*Y[jn13]*rate_eval.n_N13__p_C13 -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__n_N13 )
        )
 
+    dYdt[jc14] = (
+          -Y[jc14]*rate_eval.C14__N14__weak__wc12  +
+          ( +rho*Y[jn]*Y[jc13]*rate_eval.n_C13__C14 -Y[jc14]*rate_eval.C14__n_C13 ) +
+          ( -rho*Y[jp]*Y[jc14]*rate_eval.p_C14__N15 +Y[jn15]*rate_eval.N15__p_C14 ) +
+          ( -rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__O18 +Y[jo18]*rate_eval.O18__He4_C14 ) +
+          ( -rho*Y[jd]*Y[jc14]*rate_eval.d_C14__n_N15 +rho*Y[jn]*Y[jn15]*rate_eval.n_N15__d_C14 ) +
+          ( +rho*Y[jn]*Y[jn14]*rate_eval.n_N14__p_C14 -rho*Y[jp]*Y[jc14]*rate_eval.p_C14__n_N14 ) +
+          ( +rho*Y[jn]*Y[jo17]*rate_eval.n_O17__He4_C14 -rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__n_O17 )
+       )
+
     dYdt[jn13] = (
           -Y[jn13]*rate_eval.N13__C13__weak__wc12  +
           ( +rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13 -Y[jn13]*rate_eval.N13__p_C12 ) +
@@ -9468,6 +10064,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
        )
 
     dYdt[jn14] = (
+          +Y[jc14]*rate_eval.C14__N14__weak__wc12  +
           +Y[jo14]*rate_eval.O14__N14__weak__wc12  +
           ( +rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14 -Y[jn14]*rate_eval.N14__p_C13 ) +
           ( +rho*Y[jn]*Y[jn13]*rate_eval.n_N13__N14 -Y[jn14]*rate_eval.N14__n_N13 ) +
@@ -9476,6 +10073,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__F18 +Y[jf18]*rate_eval.F18__He4_N14 ) +
           ( +rho*Y[jhe4]*Y[jc11]*rate_eval.He4_C11__p_N14 -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__He4_C11 ) +
           ( +rho*Y[jd]*Y[jc13]*rate_eval.d_C13__n_N14 -rho*Y[jn]*Y[jn14]*rate_eval.n_N14__d_C13 ) +
+          ( -rho*Y[jn]*Y[jn14]*rate_eval.n_N14__p_C14 +rho*Y[jp]*Y[jc14]*rate_eval.p_C14__n_N14 ) +
           ( +rho*Y[jn]*Y[jo14]*rate_eval.n_O14__p_N14 -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__n_O14 ) +
           ( +rho*Y[jp]*Y[jo17]*rate_eval.p_O17__He4_N14 -rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__p_O17 ) +
           ( +rho*Y[jn]*Y[jf17]*rate_eval.n_F17__He4_N14 -rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__n_F17 )
@@ -9483,9 +10081,11 @@ def rhs_eq(t, Y, rho, T, screen_func):
 
     dYdt[jn15] = (
           +Y[jo15]*rate_eval.O15__N15__weak__wc12  +
+          ( +rho*Y[jp]*Y[jc14]*rate_eval.p_C14__N15 -Y[jn15]*rate_eval.N15__p_C14 ) +
           ( +rho*Y[jn]*Y[jn14]*rate_eval.n_N14__N15 -Y[jn15]*rate_eval.N15__n_N14 ) +
           ( -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__O16 +Y[jo16]*rate_eval.O16__p_N15 ) +
           ( -rho*Y[jhe4]*Y[jn15]*rate_eval.He4_N15__F19 +Y[jf19]*rate_eval.F19__He4_N15 ) +
+          ( +rho*Y[jd]*Y[jc14]*rate_eval.d_C14__n_N15 -rho*Y[jn]*Y[jn15]*rate_eval.n_N15__d_C14 ) +
           ( -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12 +rho*Y[jhe4]*Y[jc12]*rate_eval.He4_C12__p_N15 ) +
           ( +rho*Y[jn]*Y[jo15]*rate_eval.n_O15__p_N15 -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__n_O15 ) +
           ( +rho*Y[jp]*Y[jo18]*rate_eval.p_O18__He4_N15 -rho*Y[jhe4]*Y[jn15]*rate_eval.He4_N15__p_O18 ) +
@@ -9542,6 +10142,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
           ( -rho*Y[jn]*Y[jo17]*rate_eval.n_O17__O18 +Y[jo18]*rate_eval.O18__n_O17 ) +
           ( -rho*Y[jp]*Y[jo17]*rate_eval.p_O17__F18 +Y[jf18]*rate_eval.F18__p_O17 ) +
           ( -rho*Y[jhe4]*Y[jo17]*rate_eval.He4_O17__Ne21 +Y[jne21]*rate_eval.Ne21__He4_O17 ) +
+          ( -rho*Y[jn]*Y[jo17]*rate_eval.n_O17__He4_C14 +rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__n_O17 ) +
           ( -rho*Y[jp]*Y[jo17]*rate_eval.p_O17__He4_N14 +rho*Y[jhe4]*Y[jn14]*rate_eval.He4_N14__p_O17 ) +
           ( -rho*Y[jhe4]*Y[jo17]*rate_eval.He4_O17__n_Ne20 +rho*Y[jn]*Y[jne20]*rate_eval.n_Ne20__He4_O17 ) +
           ( +rho*Y[jn]*Y[jf17]*rate_eval.n_F17__p_O17 -rho*Y[jp]*Y[jo17]*rate_eval.p_O17__n_F17 )
@@ -9549,6 +10150,7 @@ def rhs_eq(t, Y, rho, T, screen_func):
 
     dYdt[jo18] = (
           +Y[jf18]*rate_eval.F18__O18__weak__wc12  +
+          ( +rho*Y[jhe4]*Y[jc14]*rate_eval.He4_C14__O18 -Y[jo18]*rate_eval.O18__He4_C14 ) +
           ( +rho*Y[jn]*Y[jo17]*rate_eval.n_O17__O18 -Y[jo18]*rate_eval.O18__n_O17 ) +
           ( -rho*Y[jp]*Y[jo18]*rate_eval.p_O18__F19 +Y[jf19]*rate_eval.F19__p_O18 ) +
           ( -rho*Y[jp]*Y[jo18]*rate_eval.p_O18__He4_N15 +rho*Y[jhe4]*Y[jn15]*rate_eval.He4_N15__p_O18 ) +
@@ -10029,7 +10631,9 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     n__p__weak__wc12(rate_eval, tf)
     t__He3__weak__wc12(rate_eval, tf)
     He3__t__weak__electron_capture(rate_eval, tf)
+    He6__Li6__weak__wc12(rate_eval, tf)
     Be7__Li7__weak__electron_capture(rate_eval, tf)
+    C14__N14__weak__wc12(rate_eval, tf)
     N13__C13__weak__wc12(rate_eval, tf)
     O14__N14__weak__wc12(rate_eval, tf)
     O15__N15__weak__wc12(rate_eval, tf)
@@ -10068,10 +10672,12 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     C11__He4_Be7(rate_eval, tf)
     C12__n_C11(rate_eval, tf)
     C13__n_C12(rate_eval, tf)
+    C14__n_C13(rate_eval, tf)
     N13__p_C12(rate_eval, tf)
     N14__n_N13(rate_eval, tf)
     N14__p_C13(rate_eval, tf)
     N15__n_N14(rate_eval, tf)
+    N15__p_C14(rate_eval, tf)
     O14__p_N13(rate_eval, tf)
     O15__n_O14(rate_eval, tf)
     O15__p_N14(rate_eval, tf)
@@ -10080,6 +10686,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     O16__He4_C12(rate_eval, tf)
     O17__n_O16(rate_eval, tf)
     O18__n_O17(rate_eval, tf)
+    O18__He4_C14(rate_eval, tf)
     F15__p_O14(rate_eval, tf)
     F16__n_F15(rate_eval, tf)
     F16__p_O15(rate_eval, tf)
@@ -10164,7 +10771,9 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     Cl35__He4_P31(rate_eval, tf)
     Ar35__p_Cl34(rate_eval, tf)
     Ar35__He4_S31(rate_eval, tf)
+    He6__n_n_He4(rate_eval, tf)
     Li6__n_p_He4(rate_eval, tf)
+    Be9__n_He4_He4(rate_eval, tf)
     C12__He4_He4_He4(rate_eval, tf)
     n_p__d(rate_eval, tf)
     p_p__d__weak__bet_pos_(rate_eval, tf)
@@ -10186,7 +10795,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     n_C12__C13(rate_eval, tf)
     p_C12__N13(rate_eval, tf)
     He4_C12__O16(rate_eval, tf)
+    n_C13__C14(rate_eval, tf)
     p_C13__N14(rate_eval, tf)
+    p_C14__N15(rate_eval, tf)
+    He4_C14__O18(rate_eval, tf)
     n_N13__N14(rate_eval, tf)
     p_N13__O14(rate_eval, tf)
     n_N14__N15(rate_eval, tf)
@@ -10298,16 +10910,22 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     p_Li6__He4_He3(rate_eval, tf)
     d_Li6__n_Be7(rate_eval, tf)
     d_Li6__p_Li7(rate_eval, tf)
+    He4_Li6__p_Be9(rate_eval, tf)
     p_Li7__n_Be7(rate_eval, tf)
     p_Li7__d_Li6(rate_eval, tf)
     p_Li7__He4_He4(rate_eval, tf)
+    t_Li7__n_Be9(rate_eval, tf)
     n_Be7__p_Li7(rate_eval, tf)
     n_Be7__d_Li6(rate_eval, tf)
     n_Be7__He4_He4(rate_eval, tf)
+    n_Be9__t_Li7(rate_eval, tf)
+    p_Be9__He4_Li6(rate_eval, tf)
+    He4_Be9__n_C12(rate_eval, tf)
     He4_B8__p_C11(rate_eval, tf)
     p_C11__He4_B8(rate_eval, tf)
     He4_C11__n_O14(rate_eval, tf)
     He4_C11__p_N14(rate_eval, tf)
+    n_C12__He4_Be9(rate_eval, tf)
     He4_C12__n_O15(rate_eval, tf)
     He4_C12__p_N15(rate_eval, tf)
     C12_C12__n_Mg23(rate_eval, tf)
@@ -10316,14 +10934,19 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     p_C13__n_N13(rate_eval, tf)
     d_C13__n_N14(rate_eval, tf)
     He4_C13__n_O16(rate_eval, tf)
+    p_C14__n_N14(rate_eval, tf)
+    d_C14__n_N15(rate_eval, tf)
+    He4_C14__n_O17(rate_eval, tf)
     n_N13__p_C13(rate_eval, tf)
     He4_N13__n_F16(rate_eval, tf)
     He4_N13__p_O16(rate_eval, tf)
+    n_N14__p_C14(rate_eval, tf)
     n_N14__d_C13(rate_eval, tf)
     p_N14__n_O14(rate_eval, tf)
     p_N14__He4_C11(rate_eval, tf)
     He4_N14__n_F17(rate_eval, tf)
     He4_N14__p_O17(rate_eval, tf)
+    n_N15__d_C14(rate_eval, tf)
     p_N15__n_O15(rate_eval, tf)
     p_N15__He4_C12(rate_eval, tf)
     He4_N15__n_F18(rate_eval, tf)
@@ -10346,6 +10969,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     O16_O16__n_S31(rate_eval, tf)
     O16_O16__p_P31(rate_eval, tf)
     O16_O16__He4_Si28(rate_eval, tf)
+    n_O17__He4_C14(rate_eval, tf)
     p_O17__n_F17(rate_eval, tf)
     p_O17__He4_N14(rate_eval, tf)
     He4_O17__n_Ne20(rate_eval, tf)
@@ -10512,13 +11136,17 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     He3_He3__p_p_He4(rate_eval, tf)
     d_Li7__n_He4_He4(rate_eval, tf)
     d_Be7__p_He4_He4(rate_eval, tf)
+    p_Be9__d_He4_He4(rate_eval, tf)
     n_B8__p_He4_He4(rate_eval, tf)
     n_C11__He4_He4_He4(rate_eval, tf)
     t_Li7__n_n_He4_He4(rate_eval, tf)
     He3_Li7__n_p_He4_He4(rate_eval, tf)
     t_Be7__n_p_He4_He4(rate_eval, tf)
     He3_Be7__p_p_He4_He4(rate_eval, tf)
+    p_Be9__n_p_He4_He4(rate_eval, tf)
+    n_n_He4__He6(rate_eval, tf)
     n_p_He4__Li6(rate_eval, tf)
+    n_He4_He4__Be9(rate_eval, tf)
     He4_He4_He4__C12(rate_eval, tf)
     n_p_p__p_d(rate_eval, tf)
     n_n_He4__t_t(rate_eval, tf)
@@ -10527,10 +11155,12 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     n_He4_He4__d_Li7(rate_eval, tf)
     p_He4_He4__n_B8(rate_eval, tf)
     p_He4_He4__d_Be7(rate_eval, tf)
+    d_He4_He4__p_Be9(rate_eval, tf)
     He4_He4_He4__n_C11(rate_eval, tf)
     n_n_He4_He4__t_Li7(rate_eval, tf)
     n_p_He4_He4__He3_Li7(rate_eval, tf)
     n_p_He4_He4__t_Be7(rate_eval, tf)
+    n_p_He4_He4__p_Be9(rate_eval, tf)
     p_p_He4_He4__He3_Be7(rate_eval, tf)
 
     if screen_func is not None:
@@ -10607,6 +11237,16 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         rate_eval.p_C13__N14 *= scor
         rate_eval.p_C13__n_N13 *= scor
 
+        scn_fac = ScreenFactors(1, 1, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.p_C14__N15 *= scor
+        rate_eval.p_C14__n_N14 *= scor
+
+        scn_fac = ScreenFactors(2, 4, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_C14__O18 *= scor
+        rate_eval.He4_C14__n_O17 *= scor
+
         scn_fac = ScreenFactors(1, 1, 7, 13)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.p_N13__O14 *= scor
@@ -10970,6 +11610,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_He4__n_Be7 *= scor
         rate_eval.He4_He4__p_Li7 *= scor
+        rate_eval.n_He4_He4__Be9 *= scor
         rate_eval.n_He4_He4__d_Li7 *= scor
         rate_eval.n_n_He4_He4__t_Li7 *= scor
 
@@ -10978,12 +11619,31 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         rate_eval.d_Li6__n_Be7 *= scor
         rate_eval.d_Li6__p_Li7 *= scor
 
+        scn_fac = ScreenFactors(2, 4, 3, 6)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_Li6__p_Be9 *= scor
+
         scn_fac = ScreenFactors(1, 1, 3, 7)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.p_Li7__n_Be7 *= scor
         rate_eval.p_Li7__d_Li6 *= scor
         rate_eval.p_Li7__He4_He4 *= scor
 
+        scn_fac = ScreenFactors(1, 3, 3, 7)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.t_Li7__n_Be9 *= scor
+        rate_eval.t_Li7__n_n_He4_He4 *= scor
+
+        scn_fac = ScreenFactors(1, 1, 4, 9)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.p_Be9__He4_Li6 *= scor
+        rate_eval.p_Be9__d_He4_He4 *= scor
+        rate_eval.p_Be9__n_p_He4_He4 *= scor
+
+        scn_fac = ScreenFactors(2, 4, 4, 9)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.He4_Be9__n_C12 *= scor
+
         scn_fac = ScreenFactors(2, 4, 5, 8)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_B8__p_C11 *= scor
@@ -11011,6 +11671,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_C13__n_O16 *= scor
 
+        scn_fac = ScreenFactors(1, 2, 6, 14)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.d_C14__n_N15 *= scor
+
         scn_fac = ScreenFactors(2, 4, 7, 13)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He4_N13__n_F16 *= scor
@@ -11129,10 +11793,6 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.d_Be7__p_He4_He4 *= scor
 
-        scn_fac = ScreenFactors(1, 3, 3, 7)
-        scor = screen_func(plasma_state, scn_fac)
-        rate_eval.t_Li7__n_n_He4_He4 *= scor
-
         scn_fac = ScreenFactors(2, 3, 3, 7)
         scor = screen_func(plasma_state, scn_fac)
         rate_eval.He3_Li7__n_p_He4_He4 *= scor
@@ -11162,6 +11822,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
         rate_eval.p_He4_He4__d_Be7 *= scor
         rate_eval.n_p_He4_He4__He3_Li7 *= scor
         rate_eval.n_p_He4_He4__t_Be7 *= scor
+        rate_eval.n_p_He4_He4__p_Be9 *= scor
+
+        scn_fac = ScreenFactors(1, 2, 2, 4)
+        scor = screen_func(plasma_state, scn_fac)
+        rate_eval.d_He4_He4__p_Be9 *= scor
 
     jac = np.zeros((nnuc, nnuc), dtype=np.float64)
 
@@ -11173,6 +11838,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jli6]*rate_eval.n_Li6__Li7
        -rho*Y[jc11]*rate_eval.n_C11__C12
        -rho*Y[jc12]*rate_eval.n_C12__C13
+       -rho*Y[jc13]*rate_eval.n_C13__C14
        -rho*Y[jn13]*rate_eval.n_N13__N14
        -rho*Y[jn14]*rate_eval.n_N14__N15
        -rho*Y[jo14]*rate_eval.n_O14__O15
@@ -11210,13 +11876,18 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jbe7]*rate_eval.n_Be7__p_Li7
        -rho*Y[jbe7]*rate_eval.n_Be7__d_Li6
        -rho*Y[jbe7]*rate_eval.n_Be7__He4_He4
+       -rho*Y[jbe9]*rate_eval.n_Be9__t_Li7
+       -rho*Y[jc12]*rate_eval.n_C12__He4_Be9
        -rho*Y[jn13]*rate_eval.n_N13__p_C13
+       -rho*Y[jn14]*rate_eval.n_N14__p_C14
        -rho*Y[jn14]*rate_eval.n_N14__d_C13
+       -rho*Y[jn15]*rate_eval.n_N15__d_C14
        -rho*Y[jo14]*rate_eval.n_O14__p_N14
        -rho*Y[jo14]*rate_eval.n_O14__He4_C11
        -rho*Y[jo15]*rate_eval.n_O15__p_N15
        -rho*Y[jo15]*rate_eval.n_O15__He4_C12
        -rho*Y[jo16]*rate_eval.n_O16__He4_C13
+       -rho*Y[jo17]*rate_eval.n_O17__He4_C14
        -rho*Y[jf15]*rate_eval.n_F15__p_O15
        -rho*Y[jf16]*rate_eval.n_F16__p_O16
        -rho*Y[jf16]*rate_eval.n_F16__He4_N13
@@ -11270,7 +11941,9 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jar35]*rate_eval.n_Ar35__He4_S32
        -rho*Y[jb8]*rate_eval.n_B8__p_He4_He4
        -rho*Y[jc11]*rate_eval.n_C11__He4_He4_He4
+       -2*5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__He6
        -rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6
+       -5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__Be9
        -5.00000000000000e-01*rho**2*Y[jp]**2*rate_eval.n_p_p__p_d
        -2*5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__t_t
        -rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3
@@ -11278,6 +11951,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -2*2.50000000000000e-01*rho**3*2*Y[jn]*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7
        -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        )
 
     jac[jn, jp] = (
@@ -11287,9 +11961,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho**2*Y[jn]*Y[jhe4]*rate_eval.n_p_He4__t_He3
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        +rho*Y[jt]*rate_eval.p_t__n_He3
        +rho*Y[jli7]*rate_eval.p_Li7__n_Be7
        +rho*Y[jc13]*rate_eval.p_C13__n_N13
+       +rho*Y[jc14]*rate_eval.p_C14__n_N14
        +rho*Y[jn14]*rate_eval.p_N14__n_O14
        +rho*Y[jn15]*rate_eval.p_N15__n_O15
        +rho*Y[jo15]*rate_eval.p_O15__n_F15
@@ -11314,6 +11990,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[js34]*rate_eval.p_S34__n_Cl34
        +rho*Y[jcl35]*rate_eval.p_Cl35__n_Ar35
        +rho*Y[jd]*rate_eval.p_d__n_p_p
+       +rho*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4
        +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.p_He4_He4__n_B8
        )
 
@@ -11324,6 +12001,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jt]*rate_eval.d_t__n_He4
        +rho*Y[jli6]*rate_eval.d_Li6__n_Be7
        +rho*Y[jc13]*rate_eval.d_C13__n_N14
+       +rho*Y[jc14]*rate_eval.d_C14__n_N15
        +rho*Y[jp]*rate_eval.p_d__n_p_p
        +rho*Y[jli7]*rate_eval.d_Li7__n_He4_He4
        )
@@ -11333,6 +12011,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jp]*rate_eval.p_t__n_He3
        +rho*Y[jd]*rate_eval.d_t__n_He4
        +rho*Y[jhe4]*rate_eval.He4_t__n_Li6
+       +rho*Y[jli7]*rate_eval.t_Li7__n_Be9
        +2*5.00000000000000e-01*rho*2*Y[jt]*rate_eval.t_t__n_n_He4
        +rho*Y[jhe3]*rate_eval.t_He3__n_p_He4
        +2*rho*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4
@@ -11349,19 +12028,24 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jn, jhe4] = (
        -rho*Y[jn]*rate_eval.n_He4__d_t
+       -2*5.00000000000000e-01*rho**2*Y[jn]**2*rate_eval.n_n_He4__He6
        -rho**2*Y[jn]*Y[jp]*rate_eval.n_p_He4__Li6
+       -5.00000000000000e-01*rho**2*Y[jn]*2*Y[jhe4]*rate_eval.n_He4_He4__Be9
        -2*5.00000000000000e-01*rho**2*Y[jn]**2*rate_eval.n_n_He4__t_t
        -rho**2*Y[jn]*Y[jp]*rate_eval.n_p_He4__t_He3
        -5.00000000000000e-01*rho**2*Y[jn]*2*Y[jhe4]*rate_eval.n_He4_He4__d_Li7
        -2*2.50000000000000e-01*rho**3*Y[jn]**2*2*Y[jhe4]*rate_eval.n_n_He4_He4__t_Li7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__p_Be9
        +rate_eval.He4__n_He3
        +rho*Y[jt]*rate_eval.He4_t__n_Li6
        +5.00000000000000e-01*rho*2*Y[jhe4]*rate_eval.He4_He4__n_Be7
+       +rho*Y[jbe9]*rate_eval.He4_Be9__n_C12
        +rho*Y[jc11]*rate_eval.He4_C11__n_O14
        +rho*Y[jc12]*rate_eval.He4_C12__n_O15
        +rho*Y[jc13]*rate_eval.He4_C13__n_O16
+       +rho*Y[jc14]*rate_eval.He4_C14__n_O17
        +rho*Y[jn13]*rate_eval.He4_N13__n_F16
        +rho*Y[jn14]*rate_eval.He4_N14__n_F17
        +rho*Y[jn15]*rate_eval.He4_N15__n_F18
@@ -11393,6 +12077,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +1.66666666666667e-01*rho**2*3*Y[jhe4]**2*rate_eval.He4_He4_He4__n_C11
        )
 
+    jac[jn, jhe6] = (
+       +2*rate_eval.He6__n_n_He4
+       )
+
     jac[jn, jli6] = (
        -rho*Y[jn]*rate_eval.n_Li6__Li7
        -rho*Y[jn]*rate_eval.n_Li6__He4_t
@@ -11403,6 +12091,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     jac[jn, jli7] = (
        +rate_eval.Li7__n_Li6
        +rho*Y[jp]*rate_eval.p_Li7__n_Be7
+       +rho*Y[jt]*rate_eval.t_Li7__n_Be9
        +rho*Y[jd]*rate_eval.d_Li7__n_He4_He4
        +2*rho*Y[jt]*rate_eval.t_Li7__n_n_He4_He4
        +rho*Y[jhe3]*rate_eval.He3_Li7__n_p_He4_He4
@@ -11415,6 +12104,13 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jt]*rate_eval.t_Be7__n_p_He4_He4
        )
 
+    jac[jn, jbe9] = (
+       -rho*Y[jn]*rate_eval.n_Be9__t_Li7
+       +rate_eval.Be9__n_He4_He4
+       +rho*Y[jhe4]*rate_eval.He4_Be9__n_C12
+       +rho*Y[jp]*rate_eval.p_Be9__n_p_He4_He4
+       )
+
     jac[jn, jb8] = (
        -rho*Y[jn]*rate_eval.n_B8__p_He4_He4
        )
@@ -11427,6 +12123,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jn, jc12] = (
        -rho*Y[jn]*rate_eval.n_C12__C13
+       -rho*Y[jn]*rate_eval.n_C12__He4_Be9
        +rate_eval.C12__n_C11
        +rho*Y[jhe4]*rate_eval.He4_C12__n_O15
        +5.00000000000000e-01*rho*2*Y[jc12]*rate_eval.C12_C12__n_Mg23
@@ -11434,12 +12131,20 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jn, jc13] = (
+       -rho*Y[jn]*rate_eval.n_C13__C14
        +rate_eval.C13__n_C12
        +rho*Y[jp]*rate_eval.p_C13__n_N13
        +rho*Y[jd]*rate_eval.d_C13__n_N14
        +rho*Y[jhe4]*rate_eval.He4_C13__n_O16
        )
 
+    jac[jn, jc14] = (
+       +rate_eval.C14__n_C13
+       +rho*Y[jp]*rate_eval.p_C14__n_N14
+       +rho*Y[jd]*rate_eval.d_C14__n_N15
+       +rho*Y[jhe4]*rate_eval.He4_C14__n_O17
+       )
+
     jac[jn, jn13] = (
        -rho*Y[jn]*rate_eval.n_N13__N14
        -rho*Y[jn]*rate_eval.n_N13__p_C13
@@ -11448,6 +12153,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jn, jn14] = (
        -rho*Y[jn]*rate_eval.n_N14__N15
+       -rho*Y[jn]*rate_eval.n_N14__p_C14
        -rho*Y[jn]*rate_eval.n_N14__d_C13
        +rate_eval.N14__n_N13
        +rho*Y[jp]*rate_eval.p_N14__n_O14
@@ -11455,6 +12161,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jn, jn15] = (
+       -rho*Y[jn]*rate_eval.n_N15__d_C14
        +rate_eval.N15__n_N14
        +rho*Y[jp]*rate_eval.p_N15__n_O15
        +rho*Y[jhe4]*rate_eval.He4_N15__n_F18
@@ -11486,6 +12193,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jn, jo17] = (
        -rho*Y[jn]*rate_eval.n_O17__O18
+       -rho*Y[jn]*rate_eval.n_O17__He4_C14
        +rate_eval.O17__n_O16
        +rho*Y[jp]*rate_eval.p_O17__n_F17
        +rho*Y[jhe4]*rate_eval.He4_O17__n_Ne20
@@ -11752,10 +12460,12 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3
        -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        +rate_eval.n__p__weak__wc12
        +rho*Y[jhe3]*rate_eval.n_He3__p_t
        +rho*Y[jbe7]*rate_eval.n_Be7__p_Li7
        +rho*Y[jn13]*rate_eval.n_N13__p_C13
+       +rho*Y[jn14]*rate_eval.n_N14__p_C14
        +rho*Y[jo14]*rate_eval.n_O14__p_N14
        +rho*Y[jo15]*rate_eval.n_O15__p_N15
        +rho*Y[jf15]*rate_eval.n_F15__p_O15
@@ -11781,6 +12491,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jar35]*rate_eval.n_Ar35__p_Cl35
        +rho*Y[jb8]*rate_eval.n_B8__p_He4_He4
        +5.00000000000000e-01*rho**2*Y[jp]**2*rate_eval.n_p_p__p_d
+       +5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        )
 
     jac[jp, jp] = (
@@ -11794,6 +12505,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jbe7]*rate_eval.p_Be7__B8
        -rho*Y[jc12]*rate_eval.p_C12__N13
        -rho*Y[jc13]*rate_eval.p_C13__N14
+       -rho*Y[jc14]*rate_eval.p_C14__N15
        -rho*Y[jn13]*rate_eval.p_N13__O14
        -rho*Y[jn14]*rate_eval.p_N14__O15
        -rho*Y[jn15]*rate_eval.p_N15__O16
@@ -11832,8 +12544,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jli7]*rate_eval.p_Li7__n_Be7
        -rho*Y[jli7]*rate_eval.p_Li7__d_Li6
        -rho*Y[jli7]*rate_eval.p_Li7__He4_He4
+       -rho*Y[jbe9]*rate_eval.p_Be9__He4_Li6
        -rho*Y[jc11]*rate_eval.p_C11__He4_B8
        -rho*Y[jc13]*rate_eval.p_C13__n_N13
+       -rho*Y[jc14]*rate_eval.p_C14__n_N14
        -rho*Y[jn14]*rate_eval.p_N14__n_O14
        -rho*Y[jn14]*rate_eval.p_N14__He4_C11
        -rho*Y[jn15]*rate_eval.p_N15__n_O15
@@ -11896,6 +12610,8 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jcl35]*rate_eval.p_Cl35__He4_S32
        -rho*Y[jar35]*rate_eval.p_Ar35__He4_Cl32
        -rho*Y[jd]*rate_eval.p_d__n_p_p
+       -rho*Y[jbe9]*rate_eval.p_Be9__d_He4_He4
+       -rho*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4
        -rho**2*Y[jn]*Y[jhe4]*rate_eval.n_p_He4__Li6
        -2*5.00000000000000e-01*rho**2*Y[jn]*2*Y[jp]*rate_eval.n_p_p__p_d
        -rho**2*Y[jn]*Y[jhe4]*rate_eval.n_p_He4__t_He3
@@ -11904,9 +12620,12 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        -2*2.50000000000000e-01*rho**3*2*Y[jp]*Y[jhe4]**2*rate_eval.p_p_He4_He4__He3_Be7
        +2*rho*Y[jd]*rate_eval.p_d__n_p_p
+       +rho*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4
        +5.00000000000000e-01*rho**2*Y[jn]*2*Y[jp]*rate_eval.n_p_p__p_d
+       +5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        )
 
     jac[jp, jd] = (
@@ -11918,6 +12637,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jli6]*rate_eval.d_Li6__p_Li7
        +2*rho*Y[jp]*rate_eval.p_d__n_p_p
        +rho*Y[jbe7]*rate_eval.d_Be7__p_He4_He4
+       +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9
        )
 
     jac[jp, jt] = (
@@ -11949,10 +12669,12 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -5.00000000000000e-01*rho**2*Y[jp]*2*Y[jhe4]*rate_eval.p_He4_He4__d_Be7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__He3_Li7
        -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__t_Be7
+       -5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__p_Be9
        -2*2.50000000000000e-01*rho**3*Y[jp]**2*2*Y[jhe4]*rate_eval.p_p_He4_He4__He3_Be7
        +rate_eval.He4__p_t
        +rho*Y[jhe3]*rate_eval.He4_He3__p_Li6
        +5.00000000000000e-01*rho*2*Y[jhe4]*rate_eval.He4_He4__p_Li7
+       +rho*Y[jli6]*rate_eval.He4_Li6__p_Be9
        +rho*Y[jb8]*rate_eval.He4_B8__p_C11
        +rho*Y[jc11]*rate_eval.He4_C11__p_N14
        +rho*Y[jc12]*rate_eval.He4_C12__p_N15
@@ -11988,6 +12710,8 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[js31]*rate_eval.He4_S31__p_Cl34
        +rho*Y[js32]*rate_eval.He4_S32__p_Cl35
        +rho*Y[jcl32]*rate_eval.He4_Cl32__p_Ar35
+       +5.00000000000000e-01*rho**2*Y[jd]*2*Y[jhe4]*rate_eval.d_He4_He4__p_Be9
+       +5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__p_Be9
        )
 
     jac[jp, jli6] = (
@@ -11995,6 +12719,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_Li6__He4_He3
        +rate_eval.Li6__n_p_He4
        +rho*Y[jd]*rate_eval.d_Li6__p_Li7
+       +rho*Y[jhe4]*rate_eval.He4_Li6__p_Be9
        )
 
     jac[jp, jli7] = (
@@ -12013,6 +12738,13 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +2*rho*Y[jhe3]*rate_eval.He3_Be7__p_p_He4_He4
        )
 
+    jac[jp, jbe9] = (
+       -rho*Y[jp]*rate_eval.p_Be9__He4_Li6
+       -rho*Y[jp]*rate_eval.p_Be9__d_He4_He4
+       -rho*Y[jp]*rate_eval.p_Be9__n_p_He4_He4
+       +rho*Y[jp]*rate_eval.p_Be9__n_p_He4_He4
+       )
+
     jac[jp, jb8] = (
        +rate_eval.B8__p_Be7
        +rho*Y[jhe4]*rate_eval.He4_B8__p_C11
@@ -12037,6 +12769,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_C13__n_N13
        )
 
+    jac[jp, jc14] = (
+       -rho*Y[jp]*rate_eval.p_C14__N15
+       -rho*Y[jp]*rate_eval.p_C14__n_N14
+       )
+
     jac[jp, jn13] = (
        -rho*Y[jp]*rate_eval.p_N13__O14
        +rate_eval.N13__p_C12
@@ -12049,6 +12786,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_N14__n_O14
        -rho*Y[jp]*rate_eval.p_N14__He4_C11
        +rate_eval.N14__p_C13
+       +rho*Y[jn]*rate_eval.n_N14__p_C14
        +rho*Y[jhe4]*rate_eval.He4_N14__p_O17
        )
 
@@ -12056,6 +12794,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_N15__O16
        -rho*Y[jp]*rate_eval.p_N15__n_O15
        -rho*Y[jp]*rate_eval.p_N15__He4_C12
+       +rate_eval.N15__p_C14
        +rho*Y[jhe4]*rate_eval.He4_N15__p_O18
        )
 
@@ -12364,6 +13103,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jhe4]*rate_eval.n_He4__d_t
        +rho*Y[jbe7]*rate_eval.n_Be7__d_Li6
        +rho*Y[jn14]*rate_eval.n_N14__d_C13
+       +rho*Y[jn15]*rate_eval.n_N15__d_C14
        +5.00000000000000e-01*rho**2*Y[jp]**2*rate_eval.n_p_p__p_d
        +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7
        )
@@ -12377,6 +13117,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +2*rho*Y[jt]*rate_eval.p_t__d_d
        +rho*Y[jhe4]*rate_eval.p_He4__d_He3
        +rho*Y[jli7]*rate_eval.p_Li7__d_Li6
+       +rho*Y[jbe9]*rate_eval.p_Be9__d_He4_He4
        +5.00000000000000e-01*rho**2*Y[jn]*2*Y[jp]*rate_eval.n_p_p__p_d
        +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7
        )
@@ -12395,9 +13136,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jli6]*rate_eval.d_Li6__n_Be7
        -rho*Y[jli6]*rate_eval.d_Li6__p_Li7
        -rho*Y[jc13]*rate_eval.d_C13__n_N14
+       -rho*Y[jc14]*rate_eval.d_C14__n_N15
        -rho*Y[jp]*rate_eval.p_d__n_p_p
        -rho*Y[jli7]*rate_eval.d_Li7__n_He4_He4
        -rho*Y[jbe7]*rate_eval.d_Be7__p_He4_He4
+       -5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9
        )
 
     jac[jd, jt] = (
@@ -12417,6 +13160,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     jac[jd, jhe4] = (
        -rho*Y[jd]*rate_eval.He4_d__Li6
        -rho*Y[jd]*rate_eval.d_He4__t_He3
+       -5.00000000000000e-01*rho**2*Y[jd]*2*Y[jhe4]*rate_eval.d_He4_He4__p_Be9
        +2*rate_eval.He4__d_d
        +rho*Y[jn]*rate_eval.n_He4__d_t
        +rho*Y[jp]*rate_eval.p_He4__d_He3
@@ -12440,19 +13184,32 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jn]*rate_eval.n_Be7__d_Li6
        )
 
+    jac[jd, jbe9] = (
+       +rho*Y[jp]*rate_eval.p_Be9__d_He4_He4
+       )
+
     jac[jd, jc13] = (
        -rho*Y[jd]*rate_eval.d_C13__n_N14
        )
 
+    jac[jd, jc14] = (
+       -rho*Y[jd]*rate_eval.d_C14__n_N15
+       )
+
     jac[jd, jn14] = (
        +rho*Y[jn]*rate_eval.n_N14__d_C13
        )
 
+    jac[jd, jn15] = (
+       +rho*Y[jn]*rate_eval.n_N15__d_C14
+       )
+
     jac[jt, jn] = (
        +rho*Y[jd]*rate_eval.n_d__t
        +rho*Y[jhe3]*rate_eval.n_He3__p_t
        +rho*Y[jhe4]*rate_eval.n_He4__d_t
        +rho*Y[jli6]*rate_eval.n_Li6__He4_t
+       +rho*Y[jbe9]*rate_eval.n_Be9__t_Li7
        +2*5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__t_t
        +rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3
        +2.50000000000000e-01*rho**3*2*Y[jn]*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7
@@ -12484,6 +13241,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jd]*rate_eval.d_t__n_He4
        -rho*Y[jhe4]*rate_eval.He4_t__n_Li6
        -rho*Y[jhe3]*rate_eval.t_He3__d_He4
+       -rho*Y[jli7]*rate_eval.t_Li7__n_Be9
        -2*5.00000000000000e-01*rho*2*Y[jt]*rate_eval.t_t__n_n_He4
        -rho*Y[jhe3]*rate_eval.t_He3__n_p_He4
        -rho*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4
@@ -12514,6 +13272,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jt, jli7] = (
+       -rho*Y[jt]*rate_eval.t_Li7__n_Be9
        -rho*Y[jt]*rate_eval.t_Li7__n_n_He4_He4
        +rate_eval.Li7__He4_t
        )
@@ -12522,6 +13281,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jt]*rate_eval.t_Be7__n_p_He4_He4
        )
 
+    jac[jt, jbe9] = (
+       +rho*Y[jn]*rate_eval.n_Be9__t_Li7
+       )
+
     jac[jhe3, jn] = (
        -rho*Y[jhe3]*rate_eval.n_He3__He4
        -rho*Y[jhe3]*rate_eval.n_He3__p_t
@@ -12600,19 +13363,24 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jhe4, jn] = (
        -rho*Y[jhe4]*rate_eval.n_He4__d_t
+       -5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__He6
        -rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__Li6
+       -2*5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__Be9
        -5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__t_t
        -rho**2*Y[jp]*Y[jhe4]*rate_eval.n_p_He4__t_He3
        -2*5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7
        -2*2.50000000000000e-01*rho**3*2*Y[jn]*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7
        -2*5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -2*5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -2*5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        +rho*Y[jhe3]*rate_eval.n_He3__He4
        +rho*Y[jli6]*rate_eval.n_Li6__He4_t
        +2*rho*Y[jbe7]*rate_eval.n_Be7__He4_He4
+       +rho*Y[jc12]*rate_eval.n_C12__He4_Be9
        +rho*Y[jo14]*rate_eval.n_O14__He4_C11
        +rho*Y[jo15]*rate_eval.n_O15__He4_C12
        +rho*Y[jo16]*rate_eval.n_O16__He4_C13
+       +rho*Y[jo17]*rate_eval.n_O17__He4_C14
        +rho*Y[jf16]*rate_eval.n_F16__He4_N13
        +rho*Y[jf17]*rate_eval.n_F17__He4_N14
        +rho*Y[jf18]*rate_eval.n_F18__He4_N15
@@ -12653,11 +13421,13 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -2*5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.p_He4_He4__d_Be7
        -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
        -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__t_Be7
+       -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
        -2*2.50000000000000e-01*rho**3*2*Y[jp]*Y[jhe4]**2*rate_eval.p_p_He4_He4__He3_Be7
        +rho*Y[jt]*rate_eval.p_t__He4
        +rho*Y[jhe3]*rate_eval.p_He3__He4__weak__bet_pos_
        +rho*Y[jli6]*rate_eval.p_Li6__He4_He3
        +2*rho*Y[jli7]*rate_eval.p_Li7__He4_He4
+       +rho*Y[jbe9]*rate_eval.p_Be9__He4_Li6
        +rho*Y[jc11]*rate_eval.p_C11__He4_B8
        +rho*Y[jn14]*rate_eval.p_N14__He4_C11
        +rho*Y[jn15]*rate_eval.p_N15__He4_C12
@@ -12693,11 +13463,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jcl34]*rate_eval.p_Cl34__He4_S31
        +rho*Y[jcl35]*rate_eval.p_Cl35__He4_S32
        +rho*Y[jar35]*rate_eval.p_Ar35__He4_Cl32
+       +2*rho*Y[jbe9]*rate_eval.p_Be9__d_He4_He4
+       +2*rho*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4
        )
 
     jac[jhe4, jd] = (
        -rho*Y[jhe4]*rate_eval.He4_d__Li6
        -rho*Y[jhe4]*rate_eval.d_He4__t_He3
+       -2*5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9
        +5.00000000000000e-01*rho*2*Y[jd]*rate_eval.d_d__He4
        +rho*Y[jt]*rate_eval.d_t__n_He4
        +rho*Y[jhe3]*rate_eval.d_He3__p_He4
@@ -12739,6 +13512,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jhe3]*rate_eval.He4_He3__Be7
        -rho*Y[jbe7]*rate_eval.He4_Be7__C11
        -rho*Y[jc12]*rate_eval.He4_C12__O16
+       -rho*Y[jc14]*rate_eval.He4_C14__O18
        -rho*Y[jn14]*rate_eval.He4_N14__F18
        -rho*Y[jn15]*rate_eval.He4_N15__F19
        -rho*Y[jo14]*rate_eval.He4_O14__Ne18
@@ -12776,12 +13550,15 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jd]*rate_eval.d_He4__t_He3
        -2*5.00000000000000e-01*rho*2*Y[jhe4]*rate_eval.He4_He4__n_Be7
        -2*5.00000000000000e-01*rho*2*Y[jhe4]*rate_eval.He4_He4__p_Li7
+       -rho*Y[jli6]*rate_eval.He4_Li6__p_Be9
+       -rho*Y[jbe9]*rate_eval.He4_Be9__n_C12
        -rho*Y[jb8]*rate_eval.He4_B8__p_C11
        -rho*Y[jc11]*rate_eval.He4_C11__n_O14
        -rho*Y[jc11]*rate_eval.He4_C11__p_N14
        -rho*Y[jc12]*rate_eval.He4_C12__n_O15
        -rho*Y[jc12]*rate_eval.He4_C12__p_N15
        -rho*Y[jc13]*rate_eval.He4_C13__n_O16
+       -rho*Y[jc14]*rate_eval.He4_C14__n_O17
        -rho*Y[jn13]*rate_eval.He4_N13__n_F16
        -rho*Y[jn13]*rate_eval.He4_N13__p_O16
        -rho*Y[jn14]*rate_eval.He4_N14__n_F17
@@ -12845,7 +13622,9 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[js32]*rate_eval.He4_S32__n_Ar35
        -rho*Y[js32]*rate_eval.He4_S32__p_Cl35
        -rho*Y[jcl32]*rate_eval.He4_Cl32__p_Ar35
+       -5.00000000000000e-01*rho**2*Y[jn]**2*rate_eval.n_n_He4__He6
        -rho**2*Y[jn]*Y[jp]*rate_eval.n_p_He4__Li6
+       -2*5.00000000000000e-01*rho**2*Y[jn]*2*Y[jhe4]*rate_eval.n_He4_He4__Be9
        -3*1.66666666666667e-01*rho**2*3*Y[jhe4]**2*rate_eval.He4_He4_He4__C12
        -5.00000000000000e-01*rho**2*Y[jn]**2*rate_eval.n_n_He4__t_t
        -rho**2*Y[jn]*Y[jp]*rate_eval.n_p_He4__t_He3
@@ -12853,14 +13632,21 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -2*5.00000000000000e-01*rho**2*Y[jn]*2*Y[jhe4]*rate_eval.n_He4_He4__d_Li7
        -2*5.00000000000000e-01*rho**2*Y[jp]*2*Y[jhe4]*rate_eval.p_He4_He4__n_B8
        -2*5.00000000000000e-01*rho**2*Y[jp]*2*Y[jhe4]*rate_eval.p_He4_He4__d_Be7
+       -2*5.00000000000000e-01*rho**2*Y[jd]*2*Y[jhe4]*rate_eval.d_He4_He4__p_Be9
        -3*1.66666666666667e-01*rho**2*3*Y[jhe4]**2*rate_eval.He4_He4_He4__n_C11
        -2*2.50000000000000e-01*rho**3*Y[jn]**2*2*Y[jhe4]*rate_eval.n_n_He4_He4__t_Li7
        -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__He3_Li7
        -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__t_Be7
+       -2*5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__p_Be9
        -2*2.50000000000000e-01*rho**3*Y[jp]**2*2*Y[jhe4]*rate_eval.p_p_He4_He4__He3_Be7
        )
 
+    jac[jhe4, jhe6] = (
+       +rate_eval.He6__n_n_He4
+       )
+
     jac[jhe4, jli6] = (
+       -rho*Y[jhe4]*rate_eval.He4_Li6__p_Be9
        +rate_eval.Li6__He4_d
        +rate_eval.Li6__n_p_He4
        +rho*Y[jn]*rate_eval.n_Li6__He4_t
@@ -12884,6 +13670,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +2*rho*Y[jhe3]*rate_eval.He3_Be7__p_p_He4_He4
        )
 
+    jac[jhe4, jbe9] = (
+       -rho*Y[jhe4]*rate_eval.He4_Be9__n_C12
+       +2*rate_eval.Be9__n_He4_He4
+       +rho*Y[jp]*rate_eval.p_Be9__He4_Li6
+       +2*rho*Y[jp]*rate_eval.p_Be9__d_He4_He4
+       +2*rho*Y[jp]*rate_eval.p_Be9__n_p_He4_He4
+       )
+
     jac[jhe4, jb8] = (
        -rho*Y[jhe4]*rate_eval.He4_B8__p_C11
        +2*rate_eval.B8__He4_He4__weak__wc12
@@ -12903,6 +13697,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jhe4]*rate_eval.He4_C12__n_O15
        -rho*Y[jhe4]*rate_eval.He4_C12__p_N15
        +3*rate_eval.C12__He4_He4_He4
+       +rho*Y[jn]*rate_eval.n_C12__He4_Be9
        +5.00000000000000e-01*rho*2*Y[jc12]*rate_eval.C12_C12__He4_Ne20
        +rho*Y[jo16]*rate_eval.C12_O16__He4_Mg24
        +rho*Y[jne20]*rate_eval.C12_Ne20__He4_Si28
@@ -12912,6 +13707,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jhe4]*rate_eval.He4_C13__n_O16
        )
 
+    jac[jhe4, jc14] = (
+       -rho*Y[jhe4]*rate_eval.He4_C14__O18
+       -rho*Y[jhe4]*rate_eval.He4_C14__n_O17
+       )
+
     jac[jhe4, jn13] = (
        -rho*Y[jhe4]*rate_eval.He4_N13__n_F16
        -rho*Y[jhe4]*rate_eval.He4_N13__p_O16
@@ -12958,11 +13758,13 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     jac[jhe4, jo17] = (
        -rho*Y[jhe4]*rate_eval.He4_O17__Ne21
        -rho*Y[jhe4]*rate_eval.He4_O17__n_Ne20
+       +rho*Y[jn]*rate_eval.n_O17__He4_C14
        +rho*Y[jp]*rate_eval.p_O17__He4_N14
        )
 
     jac[jhe4, jo18] = (
        -rho*Y[jhe4]*rate_eval.He4_O18__n_Ne21
+       +rate_eval.O18__He4_C14
        +rho*Y[jp]*rate_eval.p_O18__He4_N15
        )
 
@@ -13239,6 +14041,19 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jp]*rate_eval.p_Ar35__He4_Cl32
        )
 
+    jac[jhe6, jn] = (
+       +5.00000000000000e-01*rho**2*2*Y[jn]*Y[jhe4]*rate_eval.n_n_He4__He6
+       )
+
+    jac[jhe6, jhe4] = (
+       +5.00000000000000e-01*rho**2*Y[jn]**2*rate_eval.n_n_He4__He6
+       )
+
+    jac[jhe6, jhe6] = (
+       -rate_eval.He6__Li6__weak__wc12
+       -rate_eval.He6__n_n_He4
+       )
+
     jac[jli6, jn] = (
        -rho*Y[jli6]*rate_eval.n_Li6__Li7
        -rho*Y[jli6]*rate_eval.n_Li6__He4_t
@@ -13250,6 +14065,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jli6]*rate_eval.p_Li6__Be7
        -rho*Y[jli6]*rate_eval.p_Li6__He4_He3
        +rho*Y[jli7]*rate_eval.p_Li7__d_Li6
+       +rho*Y[jbe9]*rate_eval.p_Be9__He4_Li6
        +rho**2*Y[jn]*Y[jhe4]*rate_eval.n_p_He4__Li6
        )
 
@@ -13268,12 +14084,17 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jli6, jhe4] = (
+       -rho*Y[jli6]*rate_eval.He4_Li6__p_Be9
        +rho*Y[jd]*rate_eval.He4_d__Li6
        +rho*Y[jt]*rate_eval.He4_t__n_Li6
        +rho*Y[jhe3]*rate_eval.He4_He3__p_Li6
        +rho**2*Y[jn]*Y[jp]*rate_eval.n_p_He4__Li6
        )
 
+    jac[jli6, jhe6] = (
+       +rate_eval.He6__Li6__weak__wc12
+       )
+
     jac[jli6, jli6] = (
        -rate_eval.Li6__He4_d
        -rate_eval.Li6__n_p_He4
@@ -13283,6 +14104,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_Li6__He4_He3
        -rho*Y[jd]*rate_eval.d_Li6__n_Be7
        -rho*Y[jd]*rate_eval.d_Li6__p_Li7
+       -rho*Y[jhe4]*rate_eval.He4_Li6__p_Be9
        )
 
     jac[jli6, jli7] = (
@@ -13295,9 +14117,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jn]*rate_eval.n_Be7__d_Li6
        )
 
+    jac[jli6, jbe9] = (
+       +rho*Y[jp]*rate_eval.p_Be9__He4_Li6
+       )
+
     jac[jli7, jn] = (
        +rho*Y[jli6]*rate_eval.n_Li6__Li7
        +rho*Y[jbe7]*rate_eval.n_Be7__p_Li7
+       +rho*Y[jbe9]*rate_eval.n_Be9__t_Li7
        +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__d_Li7
        +2.50000000000000e-01*rho**3*2*Y[jn]*Y[jhe4]**2*rate_eval.n_n_He4_He4__t_Li7
        +5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__He3_Li7
@@ -13316,6 +14143,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jli7, jt] = (
+       -rho*Y[jli7]*rate_eval.t_Li7__n_Be9
        -rho*Y[jli7]*rate_eval.t_Li7__n_n_He4_He4
        +rho*Y[jhe4]*rate_eval.He4_t__Li7
        )
@@ -13343,6 +14171,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jp]*rate_eval.p_Li7__n_Be7
        -rho*Y[jp]*rate_eval.p_Li7__d_Li6
        -rho*Y[jp]*rate_eval.p_Li7__He4_He4
+       -rho*Y[jt]*rate_eval.t_Li7__n_Be9
        -rho*Y[jd]*rate_eval.d_Li7__n_He4_He4
        -rho*Y[jt]*rate_eval.t_Li7__n_n_He4_He4
        -rho*Y[jhe3]*rate_eval.He3_Li7__n_p_He4_He4
@@ -13353,6 +14182,10 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jn]*rate_eval.n_Be7__p_Li7
        )
 
+    jac[jli7, jbe9] = (
+       +rho*Y[jn]*rate_eval.n_Be9__t_Li7
+       )
+
     jac[jbe7, jn] = (
        -rho*Y[jbe7]*rate_eval.n_Be7__p_Li7
        -rho*Y[jbe7]*rate_eval.n_Be7__d_Li6
@@ -13423,6 +14256,57 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rate_eval.C11__He4_Be7
        )
 
+    jac[jbe9, jn] = (
+       -rho*Y[jbe9]*rate_eval.n_Be9__t_Li7
+       +rho*Y[jc12]*rate_eval.n_C12__He4_Be9
+       +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.n_He4_He4__Be9
+       +5.00000000000000e-01*rho**3*Y[jp]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
+       )
+
+    jac[jbe9, jp] = (
+       -rho*Y[jbe9]*rate_eval.p_Be9__He4_Li6
+       -rho*Y[jbe9]*rate_eval.p_Be9__d_He4_He4
+       -rho*Y[jbe9]*rate_eval.p_Be9__n_p_He4_He4
+       +5.00000000000000e-01*rho**3*Y[jn]*Y[jhe4]**2*rate_eval.n_p_He4_He4__p_Be9
+       )
+
+    jac[jbe9, jd] = (
+       +5.00000000000000e-01*rho**2*Y[jhe4]**2*rate_eval.d_He4_He4__p_Be9
+       )
+
+    jac[jbe9, jt] = (
+       +rho*Y[jli7]*rate_eval.t_Li7__n_Be9
+       )
+
+    jac[jbe9, jhe4] = (
+       -rho*Y[jbe9]*rate_eval.He4_Be9__n_C12
+       +rho*Y[jli6]*rate_eval.He4_Li6__p_Be9
+       +5.00000000000000e-01*rho**2*Y[jn]*2*Y[jhe4]*rate_eval.n_He4_He4__Be9
+       +5.00000000000000e-01*rho**2*Y[jd]*2*Y[jhe4]*rate_eval.d_He4_He4__p_Be9
+       +5.00000000000000e-01*rho**3*Y[jn]*Y[jp]*2*Y[jhe4]*rate_eval.n_p_He4_He4__p_Be9
+       )
+
+    jac[jbe9, jli6] = (
+       +rho*Y[jhe4]*rate_eval.He4_Li6__p_Be9
+       )
+
+    jac[jbe9, jli7] = (
+       +rho*Y[jt]*rate_eval.t_Li7__n_Be9
+       )
+
+    jac[jbe9, jbe9] = (
+       -rate_eval.Be9__n_He4_He4
+       -rho*Y[jn]*rate_eval.n_Be9__t_Li7
+       -rho*Y[jp]*rate_eval.p_Be9__He4_Li6
+       -rho*Y[jhe4]*rate_eval.He4_Be9__n_C12
+       -rho*Y[jp]*rate_eval.p_Be9__d_He4_He4
+       -rho*Y[jp]*rate_eval.p_Be9__n_p_He4_He4
+       )
+
+    jac[jbe9, jc12] = (
+       +rho*Y[jn]*rate_eval.n_C12__He4_Be9
+       )
+
     jac[jb8, jn] = (
        -rho*Y[jb8]*rate_eval.n_B8__p_He4_He4
        )
@@ -13503,6 +14387,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jc12, jn] = (
        -rho*Y[jc12]*rate_eval.n_C12__C13
+       -rho*Y[jc12]*rate_eval.n_C12__He4_Be9
        +rho*Y[jc11]*rate_eval.n_C11__C12
        +rho*Y[jo15]*rate_eval.n_O15__He4_C12
        +2*rho*Y[jmg23]*rate_eval.n_Mg23__C12_C12
@@ -13521,12 +14406,17 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jc12]*rate_eval.He4_C12__O16
        -rho*Y[jc12]*rate_eval.He4_C12__n_O15
        -rho*Y[jc12]*rate_eval.He4_C12__p_N15
+       +rho*Y[jbe9]*rate_eval.He4_Be9__n_C12
        +2*rho*Y[jne20]*rate_eval.He4_Ne20__C12_C12
        +rho*Y[jmg24]*rate_eval.He4_Mg24__C12_O16
        +rho*Y[jsi28]*rate_eval.He4_Si28__C12_Ne20
        +1.66666666666667e-01*rho**2*3*Y[jhe4]**2*rate_eval.He4_He4_He4__C12
        )
 
+    jac[jc12, jbe9] = (
+       +rho*Y[jhe4]*rate_eval.He4_Be9__n_C12
+       )
+
     jac[jc12, jc11] = (
        +rho*Y[jn]*rate_eval.n_C11__C12
        )
@@ -13537,6 +14427,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jn]*rate_eval.n_C12__C13
        -rho*Y[jp]*rate_eval.p_C12__N13
        -rho*Y[jhe4]*rate_eval.He4_C12__O16
+       -rho*Y[jn]*rate_eval.n_C12__He4_Be9
        -rho*Y[jhe4]*rate_eval.He4_C12__n_O15
        -rho*Y[jhe4]*rate_eval.He4_C12__p_N15
        -2*5.00000000000000e-01*rho*2*Y[jc12]*rate_eval.C12_C12__n_Mg23
@@ -13607,6 +14498,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jc13, jn] = (
+       -rho*Y[jc13]*rate_eval.n_C13__C14
        +rho*Y[jc12]*rate_eval.n_C12__C13
        +rho*Y[jn13]*rate_eval.n_N13__p_C13
        +rho*Y[jn14]*rate_eval.n_N14__d_C13
@@ -13632,12 +14524,17 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jc13, jc13] = (
        -rate_eval.C13__n_C12
+       -rho*Y[jn]*rate_eval.n_C13__C14
        -rho*Y[jp]*rate_eval.p_C13__N14
        -rho*Y[jp]*rate_eval.p_C13__n_N13
        -rho*Y[jd]*rate_eval.d_C13__n_N14
        -rho*Y[jhe4]*rate_eval.He4_C13__n_O16
        )
 
+    jac[jc13, jc14] = (
+       +rate_eval.C14__n_C13
+       )
+
     jac[jc13, jn13] = (
        +rate_eval.N13__C13__weak__wc12
        +rho*Y[jn]*rate_eval.n_N13__p_C13
@@ -13652,6 +14549,58 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jn]*rate_eval.n_O16__He4_C13
        )
 
+    jac[jc14, jn] = (
+       +rho*Y[jc13]*rate_eval.n_C13__C14
+       +rho*Y[jn14]*rate_eval.n_N14__p_C14
+       +rho*Y[jn15]*rate_eval.n_N15__d_C14
+       +rho*Y[jo17]*rate_eval.n_O17__He4_C14
+       )
+
+    jac[jc14, jp] = (
+       -rho*Y[jc14]*rate_eval.p_C14__N15
+       -rho*Y[jc14]*rate_eval.p_C14__n_N14
+       )
+
+    jac[jc14, jd] = (
+       -rho*Y[jc14]*rate_eval.d_C14__n_N15
+       )
+
+    jac[jc14, jhe4] = (
+       -rho*Y[jc14]*rate_eval.He4_C14__O18
+       -rho*Y[jc14]*rate_eval.He4_C14__n_O17
+       )
+
+    jac[jc14, jc13] = (
+       +rho*Y[jn]*rate_eval.n_C13__C14
+       )
+
+    jac[jc14, jc14] = (
+       -rate_eval.C14__N14__weak__wc12
+       -rate_eval.C14__n_C13
+       -rho*Y[jp]*rate_eval.p_C14__N15
+       -rho*Y[jhe4]*rate_eval.He4_C14__O18
+       -rho*Y[jp]*rate_eval.p_C14__n_N14
+       -rho*Y[jd]*rate_eval.d_C14__n_N15
+       -rho*Y[jhe4]*rate_eval.He4_C14__n_O17
+       )
+
+    jac[jc14, jn14] = (
+       +rho*Y[jn]*rate_eval.n_N14__p_C14
+       )
+
+    jac[jc14, jn15] = (
+       +rate_eval.N15__p_C14
+       +rho*Y[jn]*rate_eval.n_N15__d_C14
+       )
+
+    jac[jc14, jo17] = (
+       +rho*Y[jn]*rate_eval.n_O17__He4_C14
+       )
+
+    jac[jc14, jo18] = (
+       +rate_eval.O18__He4_C14
+       )
+
     jac[jn13, jn] = (
        -rho*Y[jn13]*rate_eval.n_N13__N14
        -rho*Y[jn13]*rate_eval.n_N13__p_C13
@@ -13706,6 +14655,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jn14, jn] = (
        -rho*Y[jn14]*rate_eval.n_N14__N15
+       -rho*Y[jn14]*rate_eval.n_N14__p_C14
        -rho*Y[jn14]*rate_eval.n_N14__d_C13
        +rho*Y[jn13]*rate_eval.n_N13__N14
        +rho*Y[jo14]*rate_eval.n_O14__p_N14
@@ -13717,6 +14667,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jn14]*rate_eval.p_N14__n_O14
        -rho*Y[jn14]*rate_eval.p_N14__He4_C11
        +rho*Y[jc13]*rate_eval.p_C13__N14
+       +rho*Y[jc14]*rate_eval.p_C14__n_N14
        +rho*Y[jo17]*rate_eval.p_O17__He4_N14
        )
 
@@ -13740,6 +14691,11 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jd]*rate_eval.d_C13__n_N14
        )
 
+    jac[jn14, jc14] = (
+       +rate_eval.C14__N14__weak__wc12
+       +rho*Y[jp]*rate_eval.p_C14__n_N14
+       )
+
     jac[jn14, jn13] = (
        +rho*Y[jn]*rate_eval.n_N13__N14
        )
@@ -13750,6 +14706,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jn]*rate_eval.n_N14__N15
        -rho*Y[jp]*rate_eval.p_N14__O15
        -rho*Y[jhe4]*rate_eval.He4_N14__F18
+       -rho*Y[jn]*rate_eval.n_N14__p_C14
        -rho*Y[jn]*rate_eval.n_N14__d_C13
        -rho*Y[jp]*rate_eval.p_N14__n_O14
        -rho*Y[jp]*rate_eval.p_N14__He4_C11
@@ -13783,6 +14740,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        )
 
     jac[jn15, jn] = (
+       -rho*Y[jn15]*rate_eval.n_N15__d_C14
        +rho*Y[jn14]*rate_eval.n_N14__N15
        +rho*Y[jo15]*rate_eval.n_O15__p_N15
        +rho*Y[jf18]*rate_eval.n_F18__He4_N15
@@ -13792,9 +14750,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jn15]*rate_eval.p_N15__O16
        -rho*Y[jn15]*rate_eval.p_N15__n_O15
        -rho*Y[jn15]*rate_eval.p_N15__He4_C12
+       +rho*Y[jc14]*rate_eval.p_C14__N15
        +rho*Y[jo18]*rate_eval.p_O18__He4_N15
        )
 
+    jac[jn15, jd] = (
+       +rho*Y[jc14]*rate_eval.d_C14__n_N15
+       )
+
     jac[jn15, jhe4] = (
        -rho*Y[jn15]*rate_eval.He4_N15__F19
        -rho*Y[jn15]*rate_eval.He4_N15__n_F18
@@ -13806,14 +14769,21 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        +rho*Y[jhe4]*rate_eval.He4_C12__p_N15
        )
 
+    jac[jn15, jc14] = (
+       +rho*Y[jp]*rate_eval.p_C14__N15
+       +rho*Y[jd]*rate_eval.d_C14__n_N15
+       )
+
     jac[jn15, jn14] = (
        +rho*Y[jn]*rate_eval.n_N14__N15
        )
 
     jac[jn15, jn15] = (
        -rate_eval.N15__n_N14
+       -rate_eval.N15__p_C14
        -rho*Y[jp]*rate_eval.p_N15__O16
        -rho*Y[jhe4]*rate_eval.He4_N15__F19
+       -rho*Y[jn]*rate_eval.n_N15__d_C14
        -rho*Y[jp]*rate_eval.p_N15__n_O15
        -rho*Y[jp]*rate_eval.p_N15__He4_C12
        -rho*Y[jhe4]*rate_eval.He4_N15__n_F18
@@ -14094,6 +15064,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jo17, jn] = (
        -rho*Y[jo17]*rate_eval.n_O17__O18
+       -rho*Y[jo17]*rate_eval.n_O17__He4_C14
        +rho*Y[jo16]*rate_eval.n_O16__O17
        +rho*Y[jf17]*rate_eval.n_F17__p_O17
        +rho*Y[jne20]*rate_eval.n_Ne20__He4_O17
@@ -14108,9 +15079,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
     jac[jo17, jhe4] = (
        -rho*Y[jo17]*rate_eval.He4_O17__Ne21
        -rho*Y[jo17]*rate_eval.He4_O17__n_Ne20
+       +rho*Y[jc14]*rate_eval.He4_C14__n_O17
        +rho*Y[jn14]*rate_eval.He4_N14__p_O17
        )
 
+    jac[jo17, jc14] = (
+       +rho*Y[jhe4]*rate_eval.He4_C14__n_O17
+       )
+
     jac[jo17, jn14] = (
        +rho*Y[jhe4]*rate_eval.He4_N14__p_O17
        )
@@ -14124,6 +15100,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
        -rho*Y[jn]*rate_eval.n_O17__O18
        -rho*Y[jp]*rate_eval.p_O17__F18
        -rho*Y[jhe4]*rate_eval.He4_O17__Ne21
+       -rho*Y[jn]*rate_eval.n_O17__He4_C14
        -rho*Y[jp]*rate_eval.p_O17__n_F17
        -rho*Y[jp]*rate_eval.p_O17__He4_N14
        -rho*Y[jhe4]*rate_eval.He4_O17__n_Ne20
@@ -14164,9 +15141,14 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jo18, jhe4] = (
        -rho*Y[jo18]*rate_eval.He4_O18__n_Ne21
+       +rho*Y[jc14]*rate_eval.He4_C14__O18
        +rho*Y[jn15]*rate_eval.He4_N15__p_O18
        )
 
+    jac[jo18, jc14] = (
+       +rho*Y[jhe4]*rate_eval.He4_C14__O18
+       )
+
     jac[jo18, jn15] = (
        +rho*Y[jhe4]*rate_eval.He4_N15__p_O18
        )
@@ -14177,6 +15159,7 @@ def jacobian_eq(t, Y, rho, T, screen_func):
 
     jac[jo18, jo18] = (
        -rate_eval.O18__n_O17
+       -rate_eval.O18__He4_C14
        -rho*Y[jp]*rate_eval.p_O18__F19
        -rho*Y[jp]*rate_eval.p_O18__n_F18
        -rho*Y[jp]*rate_eval.p_O18__He4_N15