diff --git a/CMakeLists.txt b/CMakeLists.txt
index 26d2cf9fc..1c17fb843 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,7 +5,7 @@ project(POET
   DESCRIPTION "A coupled reactive transport simulator")
 
 # specify the C++ standard
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 
 set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
diff --git a/bench/dolo/dol.pqi b/bench/dolo/dol.pqi
index a9706bcf1..30a6674d7 100644
--- a/bench/dolo/dol.pqi
+++ b/bench/dolo/dol.pqi
@@ -39,5 +39,5 @@ SOLUTION 4
     water 1
     temp 25
     Mg 0.002
-    Cl 0.002
+    Cl 0.004
 END
diff --git a/bench/dolo/dolo_inner_large.R b/bench/dolo/dolo_inner_large.R
index 802793581..155f75c85 100644
--- a/bench/dolo/dolo_inner_large.R
+++ b/bench/dolo/dolo_inner_large.R
@@ -16,8 +16,8 @@ bound_size <- 2
 
 diffusion_setup <- list(
     inner_boundaries = list(
-        "row" = c(200, 800, 800),
-        "col" = c(400, 1400, 1600),
+        "row" = c(400, 1400, 1600),
+        "col" = c(200, 800, 800),
         "sol_id" = c(3, 4, 4)
     ),
     alpha_x = 1e-6,
diff --git a/ext/iphreeqc b/ext/iphreeqc
index 4ec8b7006..749fdbc2e 160000
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@@ -1 +1 @@
-Subproject commit 4ec8b7006c215ad9fc1310505cd56d03ba4b17dd
+Subproject commit 749fdbc2e9478046bf3f270c70e5800637246712
diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp
index f48bd9a00..155b90569 100644
--- a/src/Chemistry/ChemistryModule.cpp
+++ b/src/Chemistry/ChemistryModule.cpp
@@ -1,12 +1,11 @@
 #include "ChemistryModule.hpp"
 
-#include "IPhreeqcPOET.hpp"
+#include "PhreeqcEngine.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
 #include <algorithm>
 #include <cassert>
-#include <cstddef>
 #include <cstdint>
 #include <memory>
 #include <mpi.h>
@@ -168,12 +167,16 @@ poet::ChemistryModule::ChemistryModule(
   this->n_cells = chem_params.total_grid_cells;
 
   if (!is_master) {
-    for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
-      this->phreeqc_instances[chem_params.pqc_ids[i]] =
-          std::make_unique<IPhreeqcPOET>(chem_params.database,
-                                         chem_params.pqc_scripts[i],
-                                         chem_params.pqc_sol_order, wp_size_);
+    for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) {
+      this->phreeqc_instances[pqc_id] =
+          std::make_unique<PhreeqcEngine>(pqc_config);
     }
+    // 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_);
+    // }
   }
 }
 
diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp
index 80be5a22f..7260422fb 100644
--- a/src/Chemistry/ChemistryModule.hpp
+++ b/src/Chemistry/ChemistryModule.hpp
@@ -11,7 +11,7 @@
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include "IPhreeqcPOET.hpp"
+#include "PhreeqcEngine.hpp"
 #include <array>
 #include <cstdint>
 #include <map>
@@ -378,7 +378,7 @@ protected:
 
   const InitialList::ChemistryInit params;
 
-  std::map<int, std::unique_ptr<IPhreeqcPOET>> phreeqc_instances;
+  std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
 };
 } // namespace poet
 
diff --git a/src/Chemistry/WorkerFunctions.cpp b/src/Chemistry/WorkerFunctions.cpp
index fb9b1f015..96a67898d 100644
--- a/src/Chemistry/WorkerFunctions.cpp
+++ b/src/Chemistry/WorkerFunctions.cpp
@@ -4,6 +4,7 @@
 
 #include "Chemistry/ChemistryDefs.hpp"
 
+#include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <iomanip>
@@ -280,8 +281,6 @@ void poet::ChemistryModule::WorkerReadDHTDump(
 void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
                                                  double dSimTime,
                                                  double dTimestep) {
-  // check if we actually need to start phreeqc
-  std::map<int, std::vector<std::size_t>> zone_mapping;
 
   for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
     if (work_package.mapping[wp_id] != CHEM_PQC) {
@@ -289,70 +288,24 @@ void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
     }
 
     auto curr_input = work_package.input[wp_id];
-    const auto pqc_id = curr_input[0];
+    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(curr_input.begin());
-
-    // remove NaNs from the input
     curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(),
                                     [](double d) { return std::isnan(d); }),
                      curr_input.end());
 
-    phreeqc_instance->queueCell(curr_input);
-    zone_mapping[pqc_id].push_back(wp_id);
-  }
+    phreeqc_instance->runCell(curr_input, dTimestep);
 
-  if (zone_mapping.empty()) {
-    return;
-  }
-
-  for (const auto &[pqc_id, eval_cells] : zone_mapping) {
-    if (eval_cells.empty()) {
-      continue;
-    }
-
-    this->phreeqc_instances[pqc_id]->runQueuedCells(dTimestep);
-
-    std::vector<std::vector<double>> pqc_out;
-    this->phreeqc_instances[pqc_id]->dequeueCells(pqc_out);
-
-    for (std::size_t i = 0; i < eval_cells.size(); i++) {
-      std::size_t wp_id = eval_cells[i];
-      std::size_t output_index = 0;
-      for (std::size_t j = 1; j < work_package.output[wp_id].size(); j++) {
-        if (!(std::isnan(work_package.output[wp_id][j]))) {
-          work_package.output[wp_id][j] = pqc_out[i][output_index++];
-        }
+    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++];
       }
     }
   }
-
-  // run the phreeqc instances
-  // for (const auto &[pqc_id, phreeqc_instance] : this->phreeqc_instances) {
-  //   // if (zone_mapping.find(pqc_id) == zone_mapping.end()) {
-  //   //   continue;
-  //   // }
-  //   phreeqc_instance->runQueuedCells(dTimestep);
-
-  //   // remap the output to the work_package
-  //   std::vector<std::vector<double>> pqc_out;
-  //   phreeqc_instance->dequeueCells(pqc_out);
-
-  //   std::size_t output_id = 0;
-
-  //   for (const auto &wp_id : zone_mapping[pqc_id]) {
-  //     std::size_t output_index = 0;
-  //     for (std::size_t i = 1; i < work_package.output[wp_id].size(); i++) {
-  //       if (!(std::isnan(work_package.output[wp_id][i]))) {
-  //         work_package.output[wp_id][i] = pqc_out[output_id][output_index++];
-  //       }
-  //     }
-  //     output_id++;
-  //   }
-  // }
 }
 
 void poet::ChemistryModule::WorkerPerfToMaster(int type,
diff --git a/src/Init/ChemistryInit.cpp b/src/Init/ChemistryInit.cpp
index e5f59380e..5eb0d01f4 100644
--- a/src/Init/ChemistryInit.cpp
+++ b/src/Init/ChemistryInit.cpp
@@ -1,12 +1,16 @@
 #include "InitialList.hpp"
 
 #include <Rcpp.h>
+#include <cstddef>
 #include <vector>
 
 namespace poet {
 
 void InitialList::initChemistry(const Rcpp::List &chem) {
-  this->pqc_sol_order = this->transport_names;
+  this->pqc_solutions = std::vector<std::string>(
+      this->transport_names.begin() + 3, this->transport_names.end());
+
+  this->pqc_solution_primaries = this->phreeqc->getSolutionPrimaries();
 
   if (chem.containsElementNamed("dht_species")) {
     this->dht_species = Rcpp::as<NamedVector<uint32_t>>(chem["dht_species"]);
@@ -27,6 +31,10 @@ void InitialList::initChemistry(const Rcpp::List &chem) {
       }
     }
   }
+
+  this->field_header =
+      Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
+  this->field_header.erase(this->field_header.begin());
 }
 
 InitialList::ChemistryInit InitialList::getChemistryInit() const {
@@ -34,11 +42,26 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const {
 
   chem_init.total_grid_cells = this->n_cols * this->n_rows;
 
-  chem_init.database = database;
-  chem_init.pqc_scripts = pqc_scripts;
-  chem_init.pqc_ids = pqc_ids;
+  // chem_init.field_header = this->field_header;
 
-  chem_init.pqc_sol_order = pqc_sol_order;
+  chem_init.database = database;
+  // chem_init.pqc_scripts = pqc_scripts;
+  // chem_init.pqc_ids = pqc_ids;
+
+  for (std::size_t i = 0; i < pqc_scripts.size(); i++) {
+    POETInitCell cell = {
+        pqc_solutions,
+        pqc_solution_primaries,
+        Rcpp::as<std::vector<std::string>>(pqc_exchanger[i]),
+        Rcpp::as<std::vector<std::string>>(pqc_kinetics[i]),
+        Rcpp::as<std::vector<std::string>>(pqc_equilibrium[i]),
+        Rcpp::as<std::vector<std::string>>(pqc_surface_comps[i]),
+        Rcpp::as<std::vector<std::string>>(pqc_surface_charges[i])};
+
+    chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell};
+  }
+
+  // chem_init.pqc_sol_order = pqc_solutions;
 
   chem_init.dht_species = dht_species;
   chem_init.interp_species = interp_species;
diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp
index 1b2406a3a..3eb2fe2fd 100644
--- a/src/Init/DiffusionInit.cpp
+++ b/src/Init/DiffusionInit.cpp
@@ -7,8 +7,8 @@
 #include <set>
 
 #include "DataStructures/Field.hpp"
-#include "IPhreeqcPOET.hpp"
 #include "InitialList.hpp"
+#include "PhreeqcInit.hpp"
 
 #include <Rcpp/Function.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
@@ -54,7 +54,7 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
 }
 
 static std::vector<std::string>
-extend_transport_names(std::unique_ptr<IPhreeqcPOET> &phreeqc,
+extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
                        const Rcpp::List &boundaries_list,
                        const Rcpp::List &inner_boundaries,
                        const std::vector<std::string> &old_trans_names) {
@@ -332,7 +332,7 @@ RcppListToBoundaryMap(const std::vector<std::string> &trans_names,
       }
 
       for (std::size_t i = 0; i < type.size(); i++) {
-        if (type[i] == tug::BC_TYPE_CONSTANT) {
+        if (static_cast<tug::BC_TYPE>(type[i]) == tug::BC_TYPE_CONSTANT) {
           bc.setBoundaryElementConstant(static_cast<tug::BC_SIDE>(tug_index), i,
                                         value[i]);
         }
diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp
index 4a3acaf65..39305dbb2 100644
--- a/src/Init/GridInit.cpp
+++ b/src/Init/GridInit.cpp
@@ -1,6 +1,7 @@
 #include "InitialList.hpp"
 
-#include <IPhreeqcPOET.hpp>
+#include "PhreeqcInit.hpp"
+
 #include <RInside.h>
 #include <Rcpp/Function.h>
 #include <Rcpp/vector/instantiation.h>
@@ -15,8 +16,8 @@
 namespace poet {
 
 static Rcpp::NumericMatrix
-pqcScriptToGrid(std::unique_ptr<IPhreeqcPOET> &phreeqc, RInside &R) {
-  IPhreeqcPOET::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
+pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) {
+  PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
 
   // add "id" to the front of the column names
 
@@ -59,9 +60,9 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
   return raws;
 }
 
-static inline uint32_t getSolutionCount(std::unique_ptr<IPhreeqcPOET> &phreeqc,
+static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
                                         const Rcpp::List &initial_grid) {
-  IPhreeqcPOET::ModulesArray mod_array;
+  PhreeqcInit::ModulesArray mod_array;
   Rcpp::Function unique_R("unique");
 
   std::vector<int> row_ids =
@@ -155,7 +156,7 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
     throw std::runtime_error("Grid size must be positive.");
   }
 
-  this->phreeqc = std::make_unique<IPhreeqcPOET>(database, script);
+  this->phreeqc = std::make_unique<PhreeqcInit>(database, script);
 
   this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
   this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
@@ -178,6 +179,11 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
 
   for (const auto &id : this->pqc_ids) {
     this->pqc_scripts.push_back(pqc_raw_dumps[id]);
+    this->pqc_exchanger.push_back(phreeqc->getExchanger(id));
+    this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id));
+    this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id));
+    this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id));
+    this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id));
   }
 }
 
diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp
index 9909906cc..e4bae7093 100644
--- a/src/Init/InitialList.cpp
+++ b/src/Init/InitialList.cpp
@@ -53,12 +53,26 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
   }
   this->database =
       Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
+  this->field_header = Rcpp::as<std::vector<std::string>>(
+      setup[static_cast<int>(ExportList::CHEM_FIELD_HEADER)]);
   this->pqc_scripts = Rcpp::as<std::vector<std::string>>(
       setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)]);
   this->pqc_ids = Rcpp::as<std::vector<int>>(
       setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]);
-  this->pqc_sol_order = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOL_ORDER)]);
+  this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
+      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
+  this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
+      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
+  this->pqc_exchanger =
+      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
+  this->pqc_kinetics =
+      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
+  this->pqc_equilibrium =
+      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
+  this->pqc_surface_comps =
+      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
+  this->pqc_surface_charges =
+      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
 
   this->dht_species = NamedVector<uint32_t>(
       setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]);
@@ -92,11 +106,25 @@ Rcpp::List InitialList::exportList() {
   out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y;
 
   out[static_cast<int>(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database);
+  out[static_cast<int>(ExportList::CHEM_FIELD_HEADER)] =
+      Rcpp::wrap(this->field_header);
   out[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)] =
       Rcpp::wrap(this->pqc_scripts);
   out[static_cast<int>(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids);
-  out[static_cast<int>(ExportList::CHEM_PQC_SOL_ORDER)] =
-      Rcpp::wrap(this->pqc_sol_order);
+  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
+      Rcpp::wrap(this->pqc_solutions);
+  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
+      Rcpp::wrap(this->pqc_solution_primaries);
+  out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
+      Rcpp::wrap(this->pqc_exchanger);
+  out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
+      Rcpp::wrap(this->pqc_kinetics);
+  out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
+      Rcpp::wrap(this->pqc_equilibrium);
+  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
+      Rcpp::wrap(this->pqc_surface_comps);
+  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
+      Rcpp::wrap(this->pqc_surface_charges);
   out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species;
   out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
       Rcpp::wrap(this->interp_species);
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index 532de6b74..bc9555329 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -2,8 +2,8 @@
 
 #include "Base/RInsidePOET.hpp"
 #include "DataStructures/NamedVector.hpp"
+#include "POETInit.hpp"
 #include <Rcpp/vector/instantiation.h>
-#include <memory>
 #include <set>
 #include <tug/Boundary.hpp>
 
@@ -13,8 +13,9 @@
 #include <cstddef>
 #include <cstdint>
 
-#include <IPhreeqcPOET.hpp>
+#include <PhreeqcInit.hpp>
 #include <map>
+#include <memory>
 #include <string>
 #include <utility>
 #include <vector>
@@ -52,9 +53,16 @@ private:
     DIFFU_ALPHA_X,
     DIFFU_ALPHA_Y,
     CHEM_DATABASE,
+    CHEM_FIELD_HEADER,
     CHEM_PQC_SCRIPTS,
     CHEM_PQC_IDS,
-    CHEM_PQC_SOL_ORDER,
+    CHEM_PQC_SOLUTIONS,
+    CHEM_PQC_SOLUTION_PRIMARY,
+    CHEM_PQC_EXCHANGER,
+    CHEM_PQC_KINETICS,
+    CHEM_PQC_EQUILIBRIUM,
+    CHEM_PQC_SURFACE_COMPS,
+    CHEM_PQC_SURFACE_CHARGES,
     CHEM_DHT_SPECIES,
     CHEM_INTERP_SPECIES,
     CHEM_HOOKS,
@@ -88,7 +96,7 @@ private:
     return GridMembersString[static_cast<std::size_t>(member)];
   }
 
-  std::unique_ptr<IPhreeqcPOET> phreeqc;
+  std::unique_ptr<PhreeqcInit> phreeqc;
 
   void prepareGrid(const Rcpp::List &grid_input);
 
@@ -178,11 +186,19 @@ private:
 
   void initChemistry(const Rcpp::List &chem_input);
 
+  std::vector<std::string> field_header;
+
   std::string database;
   std::vector<std::string> pqc_scripts;
   std::vector<int> pqc_ids;
 
-  std::vector<std::string> pqc_sol_order;
+  std::vector<std::string> pqc_solutions;
+  std::vector<std::string> pqc_solution_primaries;
+  Rcpp::List pqc_exchanger;
+  Rcpp::List pqc_kinetics;
+  Rcpp::List pqc_equilibrium;
+  Rcpp::List pqc_surface_comps;
+  Rcpp::List pqc_surface_charges;
 
   NamedVector<std::uint32_t> dht_species;
 
@@ -204,10 +220,15 @@ public:
   struct ChemistryInit {
     uint32_t total_grid_cells;
 
+    // std::vector<std::string> field_header;
+
     std::string database;
-    std::vector<std::string> pqc_scripts;
-    std::vector<int> pqc_ids;
-    std::vector<std::string> pqc_sol_order;
+    // std::vector<std::string> pqc_scripts;
+    // std::vector<int> pqc_ids;
+
+    std::map<int, POETConfig> pqc_config;
+
+    // std::vector<std::string> pqc_sol_order;
 
     NamedVector<std::uint32_t> dht_species;
     NamedVector<std::uint32_t> interp_species;
diff --git a/src/poet.cpp b/src/poet.cpp
index 7dac258af..48142cbf8 100644
--- a/src/poet.cpp
+++ b/src/poet.cpp
@@ -65,15 +65,16 @@ static void init_global_functions(RInside &R) {
 
 // HACK: this is a step back as the order and also the count of fields is
 // predefined, but it will change in the future
-void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
+// static inline void writeFieldsToR(RInside &R, const Field &trans,
+//                                   const Field &chem) {
 
-  Rcpp::DataFrame t_field(trans.asSEXP());
-  R["TMP"] = t_field;
-  R.parseEval("mysetup$state_T <- TMP");
+//   Rcpp::DataFrame t_field(trans.asSEXP());
+//   R["TMP"] = t_field;
+//   R.parseEval("mysetup$state_T <- TMP");
 
-  R["TMP"] = chem.asSEXP();
-  R.parseEval("mysetup$state_C <- TMP");
-}
+//   R["TMP"] = chem.asSEXP();
+//   R.parseEval("mysetup$state_C <- TMP");
+// }
 
 enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
 
@@ -222,7 +223,26 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
   return ParseRet::PARSER_OK;
 }
 
-static Rcpp::List RunMasterLoop(const RuntimeParameters &params,
+// HACK: this is a step back as the order and also the count of fields is
+// predefined, but it will change in the future
+void call_master_iter_end(RInside &R, const Field &trans, const Field &chem) {
+  R["TMP"] = Rcpp::wrap(trans.AsVector());
+  R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
+  R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" +
+                          std::to_string(trans.GetRequestedVecSize()) +
+                          ")), TMP_PROPS)"));
+
+  R["TMP"] = Rcpp::wrap(chem.AsVector());
+  R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
+  R.parseEval(std::string("state_C <- setNames(data.frame(matrix(TMP, nrow=" +
+                          std::to_string(chem.GetRequestedVecSize()) +
+                          ")), TMP_PROPS)"));
+  R["setup"] = *global_rt_setup;
+  R.parseEval("setup <- master_iteration_end(setup, state_T, state_C)");
+  *global_rt_setup = R["setup"];
+}
+
+static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
                                 DiffusionModule &diffusion,
                                 ChemistryModule &chem) {
 
@@ -263,13 +283,7 @@ static Rcpp::List RunMasterLoop(const RuntimeParameters &params,
     // state_C after every iteration if the cmdline option
     // --ignore-results is not given (and thus the R variable
     // store_result is TRUE)
-    {
-      const auto &trans_field = diffusion.getField().asSEXP();
-      const auto &chem_field = chem.getField().asSEXP();
-
-      *global_rt_setup = master_iteration_end_R.value()(
-          *global_rt_setup, trans_field, chem_field);
-    }
+    call_master_iter_end(R, diffusion.getField(), chem.getField());
 
     diffusion.getField().update(chem.getField());
 
@@ -400,7 +414,7 @@ int main(int argc, char *argv[]) {
 
       chemistry.masterSetField(init_list.getInitialGrid());
 
-      Rcpp::List profiling = RunMasterLoop(run_params, diffusion, chemistry);
+      Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry);
 
       MSG("finished simulation loop");
 
diff --git a/util/data_evaluation/jlFun_Eval.jl b/util/data_evaluation/jlFun_Eval.jl
new file mode 100644
index 000000000..4c95d5ced
--- /dev/null
+++ b/util/data_evaluation/jlFun_Eval.jl
@@ -0,0 +1,31 @@
+using RData, DataFrames, Plots
+
+# Load all .rds files in given directory 
+function load_data(directory)
+    files = readdir(directory)
+    # data is dictionary with iteration number as key
+    data = Dict{Int,Any}()
+    for file in files
+        if (endswith(file, ".rds") && startswith(file, "iter"))
+            # extract iteration number (iter_XXX.rds)
+            iter = parse(Int, split(split(file, "_")[2], ".")[1])
+            data[iter] = load(joinpath(directory, file))
+        end
+    end
+    return data
+end
+
+function spec_to_mat(in_df::DataFrame, spec::Symbol, cols::Int)
+    specvec = in_df[!, spec]
+    specmat = transpose(reshape(specvec, cols, :))
+
+    return specmat
+end
+
+function plot_field(data::AbstractMatrix, log::Bool=false)
+    if log
+        heatmap(1:size(data, 2), 1:size(data, 1), log10.(data), c=:viridis)
+    else
+        heatmap(1:size(data, 2), 1:size(data, 1), data, c=:viridis)
+    end
+end
\ No newline at end of file