diff --git a/ext/iphreeqc b/ext/iphreeqc
index 61214a01a..38268b4aa 160000
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@@ -1 +1 @@
-Subproject commit 61214a01ad4cf99527f657e6a41c68282e4886c4
+Subproject commit 38268b4aad03e6ce4755315f4cd690f007fa2720
diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp
index a7462b0ec..dc2d430f1 100644
--- a/src/Chemistry/ChemistryModule.cpp
+++ b/src/Chemistry/ChemistryModule.cpp
@@ -2,6 +2,7 @@
 
 #include "PhreeqcEngine.hpp"
 #include "PhreeqcMatrix.hpp"
+#include "PhreeqcRunner.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
@@ -170,16 +171,9 @@ poet::ChemistryModule::ChemistryModule(
   if (!is_master) {
     PhreeqcMatrix pqc_mat =
         PhreeqcMatrix(chem_params.database, chem_params.pqc_script);
-    for (const auto &cell_id : chem_params.pqc_ids) {
-      this->phreeqc_instances[cell_id] =
-          std::make_unique<PhreeqcEngine>(pqc_mat, cell_id);
-    }
-    // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
-    //   this->phreeqc_instances[chem_params.pqc_ids[i]] =
-    //       std::make_unique<PhreeqcWrapper>(
-    //           chem_params.database, chem_params.pqc_scripts[i],
-    //           chem_params.pqc_sol_order, chem_params.field_header, wp_size_);
-    // }
+
+    this->pqc_runner =
+        std::make_unique<PhreeqcRunner>(pqc_mat.subset(chem_params.pqc_ids));
   }
 }
 
diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp
index c06293572..d2a9fca0c 100644
--- a/src/Chemistry/ChemistryModule.hpp
+++ b/src/Chemistry/ChemistryModule.hpp
@@ -8,10 +8,11 @@
 #include "ChemistryDefs.hpp"
 
 #include "Init/InitialList.hpp"
+#include "NameDouble.h"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include "PhreeqcEngine.hpp"
+#include "PhreeqcRunner.hpp"
 #include <array>
 #include <cstdint>
 #include <map>
@@ -390,7 +391,7 @@ protected:
 
   const InitialList::ChemistryInit params;
 
-  std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
+  std::unique_ptr<PhreeqcRunner> pqc_runner;
 };
 } // namespace poet
 
diff --git a/src/Chemistry/WorkerFunctions.cpp b/src/Chemistry/WorkerFunctions.cpp
index 335132d0a..a971ddd6c 100644
--- a/src/Chemistry/WorkerFunctions.cpp
+++ b/src/Chemistry/WorkerFunctions.cpp
@@ -48,13 +48,15 @@ void poet::ChemistryModule::WorkerLoop() {
     case CHEM_FIELD_INIT: {
       ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
       if (this->ai_surrogate_enabled) {
-        this->ai_surrogate_validity_vector.resize(this->n_cells); // resize statt reserve?
+        this->ai_surrogate_validity_vector.resize(
+            this->n_cells); // resize statt reserve?
       }
       break;
     }
     case CHEM_AI_BCAST_VALIDITY: {
       // Receive the index vector of valid ai surrogate predictions
-      MPI_Bcast(&this->ai_surrogate_validity_vector.front(), this->n_cells, MPI_INT, 0, this->group_comm);
+      MPI_Bcast(&this->ai_surrogate_validity_vector.front(), this->n_cells,
+                MPI_INT, 0, this->group_comm);
       break;
     }
     case CHEM_WORK_LOOP: {
@@ -187,7 +189,6 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
     }
   }
 
-
   phreeqc_time_start = MPI_Wtime();
 
   WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
@@ -300,28 +301,19 @@ void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
                                                  double dSimTime,
                                                  double dTimestep) {
 
+  std::vector<std::vector<double>> inout_chem = work_package.input;
+  std::vector<std::size_t> to_ignore;
+
   for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
     if (work_package.mapping[wp_id] != CHEM_PQC) {
-      continue;
+      to_ignore.push_back(wp_id);
     }
+  }
+  this->pqc_runner->run(inout_chem, dTimestep, to_ignore);
 
-    auto curr_input = work_package.input[wp_id];
-    const auto pqc_id = static_cast<int>(curr_input[0]);
-
-    auto &phreeqc_instance = this->phreeqc_instances[pqc_id];
-    work_package.output[wp_id] = work_package.input[wp_id];
-
-    curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(),
-                                    [](double d) { return std::isnan(d); }),
-                     curr_input.end());
-
-    phreeqc_instance->runCell(curr_input, dTimestep);
-
-    std::size_t output_index = 0;
-    for (std::size_t i = 0; i < work_package.output[wp_id].size(); i++) {
-      if (!std::isnan(work_package.output[wp_id][i])) {
-        work_package.output[wp_id][i] = curr_input[output_index++];
-      }
+  for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
+    if (work_package.mapping[wp_id] == CHEM_PQC) {
+      work_package.output[wp_id] = inout_chem[wp_id];
     }
   }
 }