diff --git a/ext/iphreeqc b/ext/iphreeqc
index d6e713074..f69561987 160000
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@@ -1 +1 @@
-Subproject commit d6e7130746f9823e0bc5f711179e676e27ba3737
+Subproject commit f695619875fcf43e411af67d592c05b199f41a15
diff --git a/ext/tug b/ext/tug
index 8c0687a6c..3ffa0ef62 160000
--- a/ext/tug
+++ b/ext/tug
@@ -1 +1 @@
-Subproject commit 8c0687a6cd4a10a79c7a554083a35eda11cc55f0
+Subproject commit 3ffa0ef6242d467a269be83751c857cb24ae532e
diff --git a/src/Base/SimParams.cpp b/src/Base/SimParams.cpp
index 3ae76dd65..8264b6c60 100644
--- a/src/Base/SimParams.cpp
+++ b/src/Base/SimParams.cpp
@@ -29,7 +29,6 @@
 
 #include <cstdint>
 #include <string>
-#include <vector>
 
 #include "Base/RInsidePOET.hpp"
 #include "argh.hpp" // Argument handler https://github.com/adishavit/argh
diff --git a/src/Base/SimParams.hpp b/src/Base/SimParams.hpp
index b4e20f262..74d6d2078 100644
--- a/src/Base/SimParams.hpp
+++ b/src/Base/SimParams.hpp
@@ -47,36 +47,37 @@ enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
  * @brief Defining all simulation parameters
  *
  */
-struct RuntimeParameters {
+// struct RuntimeParameters {
 
-  /** Count of processes in MPI_COMM_WORLD */
-  int world_size;
-  /** rank of proces in MPI_COMM_WORLD */
-  int world_rank;
-  /** indicates if DHT should be used */
-  bool dht_enabled;
-  /** apply logarithm to key before rounding */
-  bool dht_log;
-  /** indicates if timestep dt differs between iterations */
-  bool dt_differ;
-  /** Indicates, when a DHT snapshot should be written */
-  int dht_snaps;
-  /** <b>not implemented</b>: How a DHT is distributed over processes */
-  int dht_strategy;
-  /** Size of DHt per process in byter */
-  unsigned int dht_size_per_process;
-  /** Default significant digit for rounding */
-  int dht_significant_digits;
-  /** Default work package size */
-  unsigned int wp_size;
-  /** indicates if resulting grid should be stored after every iteration */
-  bool store_result;
-  /** indicating whether the progress bar during chemistry simulation should be
-   * printed or not */
-  bool print_progressbar;
+//   /** Count of processes in MPI_COMM_WORLD */
+//   int world_size;
+//   /** rank of proces in MPI_COMM_WORLD */
+//   int world_rank;
+//   /** indicates if DHT should be used */
+//   bool dht_enabled;
+//   /** apply logarithm to key before rounding */
+//   bool dht_log;
+//   /** indicates if timestep dt differs between iterations */
+//   bool dt_differ;
+//   /** Indicates, when a DHT snapshot should be written */
+//   int dht_snaps;
+//   /** <b>not implemented</b>: How a DHT is distributed over processes */
+//   int dht_strategy;
+//   /** Size of DHt per process in byter */
+//   unsigned int dht_size_per_process;
+//   /** Default significant digit for rounding */
+//   int dht_significant_digits;
+//   /** Default work package size */
+//   unsigned int wp_size;
+//   /** indicates if resulting grid should be stored after every iteration */
+//   bool store_result;
+//   /** indicating whether the progress bar during chemistry simulation should
+//   be
+//    * printed or not */
+//   bool print_progressbar;
 
-  bool interp_enabled;
-};
+//   bool interp_enabled;
+// };
 
 // using GridParams = struct s_GridParams {
 //   std::array<uint32_t, 2> n_cells;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 2e4430848..c4e8d759d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -8,6 +8,14 @@ target_sources(POETLib
   Init/ChemistryInit.cpp
   DataStructures/Field.cpp
   Transport/DiffusionModule.cpp
+  Chemistry/ChemistryModule.cpp
+  Chemistry/MasterFunctions.cpp
+  Chemistry/WorkerFunctions.cpp
+  Chemistry/SurrogateModels/DHT_Wrapper.cpp
+  Chemistry/SurrogateModels/DHT.c
+  Chemistry/SurrogateModels/HashFunctions.cpp
+  Chemistry/SurrogateModels/InterpolationModule.cpp
+  Chemistry/SurrogateModels/ProximityHashTable.cpp
 )
 
 target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
@@ -16,6 +24,7 @@ target_link_libraries(
   PUBLIC RRuntime
   PUBLIC IPhreeqcPOET
   PUBLIC tug
+  PUBLIC MPI::MPI_C
 )
 
 # add_library(poetlib
diff --git a/src/Chemistry/ChemistryDefs.hpp b/src/Chemistry/ChemistryDefs.hpp
new file mode 100644
index 000000000..71c82edf6
--- /dev/null
+++ b/src/Chemistry/ChemistryDefs.hpp
@@ -0,0 +1,23 @@
+#pragma once
+
+#include <cstdint>
+#include <vector>
+
+namespace poet {
+
+enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
+enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
+
+struct WorkPackage {
+  std::size_t size;
+  std::vector<std::vector<double>> input;
+  std::vector<std::vector<double>> output;
+  std::vector<std::uint8_t> mapping;
+
+  WorkPackage(std::size_t _size) : size(_size) {
+    input.resize(size);
+    output.resize(size);
+    mapping.resize(size, CHEM_PQC);
+  }
+};
+} // namespace poet
\ No newline at end of file
diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp
index 46531f2c6..8623ef38f 100644
--- a/src/Chemistry/ChemistryModule.cpp
+++ b/src/Chemistry/ChemistryModule.cpp
@@ -1,20 +1,16 @@
 #include "ChemistryModule.hpp"
 
+#include "IPhreeqcPOET.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include <PhreeqcRM.h>
-
 #include <algorithm>
 #include <cassert>
 #include <cstddef>
 #include <cstdint>
-#include <map>
 #include <memory>
 #include <mpi.h>
-#include <stdexcept>
 #include <string>
-#include <utility>
 #include <vector>
 
 constexpr uint32_t MB_FACTOR = 1E6;
@@ -158,25 +154,26 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
   return results;
 }
 
-poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size,
-                                       std::uint32_t maxiter,
-                                       const ChemistryParams &chem_param,
-                                       MPI_Comm communicator)
-    : PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size),
-      params(chem_param) {
+poet::ChemistryModule::ChemistryModule(
+    uint32_t wp_size_, const InitialList::ChemistryInit &chem_params,
+    MPI_Comm communicator)
+    : params(chem_params), wp_size(wp_size_), group_comm(communicator) {
+  MPI_Comm_rank(communicator, &comm_rank);
+  MPI_Comm_size(communicator, &comm_size);
 
-  MPI_Comm_size(communicator, &this->comm_size);
-  MPI_Comm_rank(communicator, &this->comm_rank);
+  this->is_sequential = comm_size == 1;
+  this->is_master = comm_rank == 0;
 
-  this->is_sequential = (this->comm_size == 1);
-  this->is_master = (this->comm_rank == 0);
+  this->n_cells = chem_params.total_grid_cells;
 
-  if (!is_sequential && is_master) {
-    MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm);
-    MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm);
+  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_);
+    }
   }
-
-  this->file_pad = std::ceil(std::log10(maxiter + 1));
 }
 
 poet::ChemistryModule::~ChemistryModule() {
@@ -185,206 +182,45 @@ poet::ChemistryModule::~ChemistryModule() {
   }
 }
 
-poet::ChemistryModule
-poet::ChemistryModule::createWorker(MPI_Comm communicator,
-                                    const ChemistryParams &chem_param) {
-  uint32_t wp_size;
-  MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, communicator);
-
-  std::uint32_t maxiter;
-  MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, communicator);
-
-  return ChemistryModule(wp_size, wp_size, maxiter, chem_param, communicator);
-}
-
-void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
-  if (this->is_master) {
-    int f_type = CHEM_INIT;
-    PropagateFunctionType(f_type);
-
-    int count = input_script_path.size();
-    ChemBCast(&count, 1, MPI_INT);
-    ChemBCast(const_cast<char *>(input_script_path.data()), count, MPI_CHAR);
-  }
-
-  this->RunFile(true, true, false, input_script_path);
-  this->RunString(true, false, false, "DELETE; -all; PRINT; -warnings 0;");
-
-  this->FindComponents();
-
-  PhreeqcRM::initializePOET(this->speciesPerModule, this->prop_names);
-  this->prop_count = prop_names.size();
-
-  char exchange = (speciesPerModule[1] == 0 ? -1 : 1);
-  char kinetics = (speciesPerModule[2] == 0 ? -1 : 1);
-  char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1);
-  char surface = (speciesPerModule[4] == 0 ? -1 : 1);
-
-  std::vector<int> ic1;
-  ic1.resize(this->nxyz * 7, -1);
-  // TODO: hardcoded reaction modules
-  for (int i = 0; i < nxyz; i++) {
-    ic1[i] = 1;                   // Solution 1
-    ic1[nxyz + i] = equilibrium;  // Equilibrium 1
-    ic1[2 * nxyz + i] = exchange; // Exchange none
-    ic1[3 * nxyz + i] = surface;  // Surface none
-    ic1[4 * nxyz + i] = -1;       // Gas phase none
-    ic1[5 * nxyz + i] = -1;       // Solid solutions none
-    ic1[6 * nxyz + i] = kinetics; // Kinetics 1
-  }
-
-  this->InitialPhreeqc2Module(ic1);
-}
-
-void poet::ChemistryModule::initializeField(const Field &trans_field) {
-
-  if (is_master) {
-    int f_type = CHEM_INIT_SPECIES;
-    PropagateFunctionType(f_type);
-  }
-
-  std::vector<std::string> essentials_backup{
-      prop_names.begin() + speciesPerModule[0], prop_names.end()};
-
-  std::vector<std::string> new_solution_names{
-      this->prop_names.begin(), this->prop_names.begin() + speciesPerModule[0]};
-
-  if (is_master) {
-    for (auto &prop : trans_field.GetProps()) {
-      if (std::find(new_solution_names.begin(), new_solution_names.end(),
-                    prop) == new_solution_names.end()) {
-        int size = prop.size();
-        ChemBCast(&size, 1, MPI_INT);
-        ChemBCast(prop.data(), prop.size(), MPI_CHAR);
-        new_solution_names.push_back(prop);
-      }
-    }
-    int end = 0;
-    ChemBCast(&end, 1, MPI_INT);
-  } else {
-    constexpr int MAXSIZE = 128;
-    MPI_Status status;
-    int recv_size;
-    char recv_buffer[MAXSIZE];
-    while (1) {
-      ChemBCast(&recv_size, 1, MPI_INT);
-      if (recv_size == 0) {
-        break;
-      }
-      ChemBCast(recv_buffer, recv_size, MPI_CHAR);
-      recv_buffer[recv_size] = '\0';
-      new_solution_names.push_back(std::string(recv_buffer));
-    }
-  }
-
-  // now sort the new values
-  std::sort(new_solution_names.begin() + 3, new_solution_names.end());
-  this->SetPOETSolutionNames(new_solution_names);
-  this->speciesPerModule[0] = new_solution_names.size();
-
-  // and append other processes than solutions
-  std::vector<std::string> new_prop_names = new_solution_names;
-  new_prop_names.insert(new_prop_names.end(), essentials_backup.begin(),
-                        essentials_backup.end());
-
-  std::vector<std::string> old_prop_names{this->prop_names};
-  this->prop_names = std::move(new_prop_names);
-  this->prop_count = prop_names.size();
-
-  if (is_master) {
-    this->n_cells = trans_field.GetRequestedVecSize();
-
-    std::vector<std::vector<double>> phreeqc_dump(this->nxyz);
-    this->getDumpedField(phreeqc_dump);
-
-    if (is_sequential) {
-      std::vector<double> init_vec;
-      for (std::size_t i = 0; i < n_cells; i++) {
-        init_vec.insert(init_vec.end(), phreeqc_dump[i].begin(),
-                        phreeqc_dump[i].end());
-      }
-
-      const auto tmp_buffer{init_vec};
-      this->unshuffleField(tmp_buffer, n_cells, prop_count, 1, init_vec);
-      this->chem_field = Field(n_cells, init_vec, prop_names);
-      return;
-    }
-
-    std::vector<double> &phreeqc_init = phreeqc_dump[0];
-    std::vector<std::vector<double>> initial_values;
-
-    for (const auto &vec : trans_field.As2DVector()) {
-      initial_values.push_back(vec);
-    }
-
-    this->base_totals = {initial_values.at(0).at(0),
-                         initial_values.at(1).at(0)};
-    ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
-
-    for (int i = speciesPerModule[0]; i < phreeqc_init.size(); i++) {
-      std::vector<double> init(n_cells, phreeqc_init[i]);
-      initial_values.push_back(std::move(init));
-    }
-
-    this->chem_field = Field(n_cells, initial_values, prop_names);
-  } else {
-    ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
-  }
-
-  if (this->params.use_dht || this->params.use_interp) {
-    initializeDHT(this->params.dht_size, this->params.dht_signifs);
-    setDHTSnapshots(this->params.dht_snaps, this->params.dht_outdir);
-    setDHTReadFile(this->params.dht_file);
-
-    this->dht_enabled = this->params.use_dht;
-
-    if (this->params.use_interp) {
-      initializeInterp(this->params.pht_max_entries, this->params.pht_size,
-                       this->params.interp_min_entries,
-                       this->params.pht_signifs);
-      this->interp_enabled = this->params.use_interp;
-    }
-  }
-}
-
 void poet::ChemistryModule::initializeDHT(
     uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
-  constexpr uint32_t MB_FACTOR = 1E6;
+  // constexpr uint32_t MB_FACTOR = 1E6;
 
-  this->dht_enabled = true;
+  // this->dht_enabled = true;
 
-  MPI_Comm dht_comm;
+  // MPI_Comm dht_comm;
 
-  if (this->is_master) {
-    MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank, &dht_comm);
-    return;
-  }
+  // if (this->is_master) {
+  //   MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank,
+  //   &dht_comm); return;
+  // }
 
-  if (!this->is_master) {
+  // if (!this->is_master) {
 
-    MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
+  //   MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
 
-    auto map_copy = key_species;
+  //   auto map_copy = key_species;
 
-    if (key_species.empty()) {
-      std::vector<std::uint32_t> default_signif(
-          this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
-      map_copy = NamedVector<std::uint32_t>(this->prop_names, default_signif);
-    }
+  //   if (key_species.empty()) {
+  //     std::vector<std::uint32_t> default_signif(
+  //         this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
+  //     map_copy = NamedVector<std::uint32_t>(this->prop_names,
+  //     default_signif);
+  //   }
 
-    auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
+  //   auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
 
-    if (this->dht) {
-      delete this->dht;
-    }
+  //   if (this->dht) {
+  //     delete this->dht;
+  //   }
 
-    const std::uint64_t dht_size = size_mb * MB_FACTOR;
+  //   const std::uint64_t dht_size = size_mb * MB_FACTOR;
 
-    this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
-                                this->prop_names, params.hooks,
-                                this->prop_count, params.use_interp);
-    this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
-  }
+  //   this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
+  //                               this->prop_names, params.hooks,
+  //                               this->prop_count, params.use_interp);
+  //   this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
+  // }
 }
 
 inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
@@ -457,42 +293,42 @@ void poet::ChemistryModule::initializeInterp(
     std::uint32_t bucket_size, std::uint32_t size_mb, std::uint32_t min_entries,
     const NamedVector<std::uint32_t> &key_species) {
 
-  if (!this->is_master) {
+  // if (!this->is_master) {
 
-    constexpr uint32_t MB_FACTOR = 1E6;
+  //   constexpr uint32_t MB_FACTOR = 1E6;
 
-    assert(this->dht);
+  //   assert(this->dht);
 
-    this->interp_enabled = true;
+  //   this->interp_enabled = true;
 
-    auto map_copy = key_species;
+  //   auto map_copy = key_species;
 
-    if (key_species.empty()) {
-      map_copy = this->dht->getKeySpecies();
-      for (std::size_t i = 0; i < map_copy.size(); i++) {
-        const std::uint32_t signif =
-            map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
-                               ? InterpolationModule::COARSE_DIFF
-                               : 0);
-        map_copy[i] = signif;
-      }
-    }
+  //   if (key_species.empty()) {
+  //     map_copy = this->dht->getKeySpecies();
+  //     for (std::size_t i = 0; i < map_copy.size(); i++) {
+  //       const std::uint32_t signif =
+  //           map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
+  //                              ? InterpolationModule::COARSE_DIFF
+  //                              : 0);
+  //       map_copy[i] = signif;
+  //     }
+  //   }
 
-    auto key_indices =
-        parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
+  //   auto key_indices =
+  //       parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
 
-    if (this->interp) {
-      this->interp.reset();
-    }
+  //   if (this->interp) {
+  //     this->interp.reset();
+  //   }
 
-    const uint64_t pht_size = size_mb * MB_FACTOR;
+  //   const uint64_t pht_size = size_mb * MB_FACTOR;
 
-    interp = std::make_unique<poet::InterpolationModule>(
-        bucket_size, pht_size, min_entries, *(this->dht), map_copy, key_indices,
-        this->prop_names, this->params.hooks);
+  //   interp = std::make_unique<poet::InterpolationModule>(
+  //       bucket_size, pht_size, min_entries, *(this->dht), map_copy,
+  //       key_indices, this->prop_names, this->params.hooks);
 
-    interp->setInterpolationFunction(inverseDistanceWeighting);
-  }
+  //   interp->setInterpolationFunction(inverseDistanceWeighting);
+  // }
 }
 
 std::vector<double>
diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp
index d343ffa31..9ca8394ea 100644
--- a/src/Chemistry/ChemistryModule.hpp
+++ b/src/Chemistry/ChemistryModule.hpp
@@ -2,16 +2,17 @@
 #ifndef CHEMISTRYMODULE_H_
 #define CHEMISTRYMODULE_H_
 
-#include "../Base/SimParams.hpp"
-#include "../DataStructures/DataStructures.hpp"
+#include "DataStructures/Field.hpp"
+#include "DataStructures/NamedVector.hpp"
+
+#include "ChemistryDefs.hpp"
+
+#include "Init/InitialList.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include <IrmResult.h>
-#include <PhreeqcRM.h>
-
+#include "IPhreeqcPOET.hpp"
 #include <array>
-#include <cstddef>
 #include <cstdint>
 #include <map>
 #include <memory>
@@ -24,7 +25,7 @@ namespace poet {
  * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
  * easy access.
  */
-class ChemistryModule : public PhreeqcRM {
+class ChemistryModule {
 public:
   /**
    * Creates a new instance of Chemistry module with given grid cell count, work
@@ -41,45 +42,40 @@ public:
    * \param wp_size Count of grid cells to fill each work package at maximum.
    * \param communicator MPI communicator to distribute work in.
    */
-  ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
-                  const ChemistryParams &chem_param, MPI_Comm communicator);
+  ChemistryModule(uint32_t wp_size,
+                  const InitialList::ChemistryInit &chem_params,
+                  MPI_Comm communicator);
 
   /**
    * Deconstructor, which frees DHT data structure if used.
    */
   ~ChemistryModule();
 
-  /**
-   * Parses the input script and extract information needed during runtime.
-   *
-   * **Only run by master**.
-   *
-   * Database must be loaded beforehand.
-   *
-   * \param input_script_path Path to input script to parse.
-   */
-  void RunInitFile(const std::string &input_script_path);
-
+  void masterSetField(Field field);
   /**
    * Run the chemical simulation with parameters set.
    */
-  void RunCells();
+  void simulate(double dt);
 
   /**
    * Returns the chemical field.
    */
-  auto GetField() const { return this->chem_field; }
+  auto &GetField() { return this->chem_field; }
   /**
    * Returns all known species names, including not only aqueous species, but
    * also equilibrium, exchange, surface and kinetic reactants.
    */
-  auto GetPropNames() const { return this->prop_names; }
+  // auto GetPropNames() const { return this->prop_names; }
 
   /**
    * Return the accumulated runtime in seconds for chemical simulation.
    */
   auto GetChemistryTime() const { return this->chem_t; }
 
+  void setFilePadding(std::uint32_t maxiter) {
+    this->file_pad = std::ceil(std::log10(maxiter + 1));
+  }
+
   /**
    * Create a new worker instance inside given MPI communicator.
    *
@@ -118,17 +114,6 @@ public:
     DHT_SNAPS_ITEREND       //!< output snapshots after each iteration
   };
 
-  /**
-   * **This function has to be run!**
-   *
-   * Merge initial values from existing module with the chemistry module and set
-   * according internal variables.
-   *
-   * \param other Field to merge chemistry with. Most likely it is something
-   * like the diffusion field.
-   */
-  void initializeField(const Field &other);
-
   /**
    * **Only called by workers!** Start the worker listening loop.
    */
@@ -243,9 +228,7 @@ protected:
                         const NamedVector<std::uint32_t> &key_species);
 
   enum {
-    CHEM_INIT,
-    CHEM_WP_SIZE,
-    CHEM_INIT_SPECIES,
+    CHEM_FIELD_INIT,
     CHEM_DHT_ENABLE,
     CHEM_DHT_SIGNIF_VEC,
     CHEM_DHT_SNAPS,
@@ -294,7 +277,7 @@ protected:
   using worker_list_t = std::vector<struct worker_info_s>;
   using workpointer_t = std::vector<double>::iterator;
 
-  void MasterRunParallel();
+  void MasterRunParallel(double dt);
   void MasterRunSequential();
 
   void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
@@ -320,8 +303,8 @@ protected:
   void WorkerPerfToMaster(int type, const struct worker_s &timings);
   void WorkerMetricsToMaster(int type);
 
-  IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
-                                  double dTimestep);
+  void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
+                            double dTimestep);
 
   std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
                                                uint32_t wp_size) const;
@@ -372,21 +355,18 @@ protected:
 
   bool print_progessbar{false};
 
-  std::uint32_t file_pad;
+  std::uint32_t file_pad{1};
 
   double chem_t{0.};
 
   uint32_t n_cells = 0;
   uint32_t prop_count = 0;
-  std::vector<std::string> prop_names;
 
   Field chem_field;
 
-  static constexpr int MODULE_COUNT = 5;
+  const InitialList::ChemistryInit &params;
 
-  const ChemistryParams &params;
-
-  std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
+  std::map<int, std::unique_ptr<IPhreeqcPOET>> phreeqc_instances;
 };
 } // namespace poet
 
diff --git a/src/Chemistry/MasterFunctions.cpp b/src/Chemistry/MasterFunctions.cpp
index 543cd21a8..9ee52be82 100644
--- a/src/Chemistry/MasterFunctions.cpp
+++ b/src/Chemistry/MasterFunctions.cpp
@@ -1,12 +1,9 @@
 #include "ChemistryModule.hpp"
 
-#include <PhreeqcRM.h>
-
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <mpi.h>
-#include <stdexcept>
 #include <vector>
 
 std::vector<uint32_t>
@@ -308,45 +305,46 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
   }
 }
 
-void poet::ChemistryModule::RunCells() {
+void poet::ChemistryModule::simulate(double dt) {
   double start_t{MPI_Wtime()};
   if (this->is_sequential) {
     MasterRunSequential();
     return;
   }
 
-  MasterRunParallel();
+  MasterRunParallel(dt);
   double end_t{MPI_Wtime()};
   this->chem_t += end_t - start_t;
 }
 
 void poet::ChemistryModule::MasterRunSequential() {
-  std::vector<double> shuffled_field =
-      shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
+  // std::vector<double> shuffled_field =
+  //     shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
 
-  std::vector<std::vector<double>> input;
-  for (std::size_t i = 0; i < n_cells; i++) {
-    input.push_back(
-        std::vector<double>(shuffled_field.begin() + (i * prop_count),
-                            shuffled_field.begin() + ((i + 1) * prop_count)));
-  }
+  // std::vector<std::vector<double>> input;
+  // for (std::size_t i = 0; i < n_cells; i++) {
+  //   input.push_back(
+  //       std::vector<double>(shuffled_field.begin() + (i * prop_count),
+  //                           shuffled_field.begin() + ((i + 1) *
+  //                           prop_count)));
+  // }
 
-  this->setDumpedField(input);
-  PhreeqcRM::RunCells();
-  this->getDumpedField(input);
+  // this->setDumpedField(input);
+  // PhreeqcRM::RunCells();
+  // this->getDumpedField(input);
 
-  shuffled_field.clear();
-  for (std::size_t i = 0; i < n_cells; i++) {
-    shuffled_field.insert(shuffled_field.end(), input[i].begin(),
-                          input[i].end());
-  }
+  // shuffled_field.clear();
+  // for (std::size_t i = 0; i < n_cells; i++) {
+  //   shuffled_field.insert(shuffled_field.end(), input[i].begin(),
+  //                         input[i].end());
+  // }
 
-  std::vector<double> out_vec{shuffled_field};
-  unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
-  chem_field = out_vec;
+  // std::vector<double> out_vec{shuffled_field};
+  // unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
+  // chem_field = out_vec;
 }
 
-void poet::ChemistryModule::MasterRunParallel() {
+void poet::ChemistryModule::MasterRunParallel(double dt) {
   /* declare most of the needed variables here */
   double seq_a, seq_b, seq_c, seq_d;
   double worker_chemistry_a, worker_chemistry_b;
@@ -360,7 +358,6 @@ void poet::ChemistryModule::MasterRunParallel() {
 
   MPI_Barrier(this->group_comm);
 
-  double dt = this->PhreeqcRM::GetTimeStep();
   static uint32_t iteration = 0;
 
   /* start time measurement of sequential part */
@@ -466,3 +463,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
 
   return wp_sizes_vector;
 }
+
+void poet::ChemistryModule::masterSetField(Field field) {
+  this->chem_field = field;
+  this->prop_count = field.GetProps().size();
+
+  int ftype = CHEM_FIELD_INIT;
+  PropagateFunctionType(ftype);
+
+  ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
+}
\ No newline at end of file
diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
index 69491b2c2..53ef72826 100644
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
+++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
@@ -22,13 +22,14 @@
 
 #include "DHT_Wrapper.hpp"
 
+#include "Rounding.hpp"
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <cstring>
-#include <iostream>
 #include <stdexcept>
 #include <vector>
 
diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
index 791be2d1b..8b741f7e1 100644
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
+++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
@@ -23,19 +23,18 @@
 #ifndef DHT_WRAPPER_H
 #define DHT_WRAPPER_H
 
-#include "../../Base/RInsidePOET.hpp"
+#include "Base/RInsidePOET.hpp"
+#include "DataStructures/NamedVector.hpp"
+
 #include "../../Base/SimParams.hpp"
-#include "../../DataStructures/DataStructures.hpp"
-#include "../enums.hpp"
-#include "HashFunctions.hpp"
+#include "Chemistry/ChemistryDefs.hpp"
+
 #include "LookupKey.hpp"
-#include "Rounding.hpp"
 
 #include <array>
 #include <cstdint>
 #include <limits>
 #include <string>
-#include <unordered_map>
 #include <utility>
 #include <vector>
 
diff --git a/src/Chemistry/SurrogateModels/Interpolation.hpp b/src/Chemistry/SurrogateModels/Interpolation.hpp
index 15361f4bf..72b0ea3bd 100644
--- a/src/Chemistry/SurrogateModels/Interpolation.hpp
+++ b/src/Chemistry/SurrogateModels/Interpolation.hpp
@@ -3,14 +3,12 @@
 #ifndef INTERPOLATION_H_
 #define INTERPOLATION_H_
 
-#include "../../Base/SimParams.hpp"
-#include "../../DataStructures/DataStructures.hpp"
+#include "DataStructures/NamedVector.hpp"
+
 #include "DHT_Wrapper.hpp"
 #include "LookupKey.hpp"
 #include "Rounding.hpp"
 
-#include <cassert>
-#include <iostream>
 #include <list>
 #include <memory>
 #include <mpi.h>
@@ -22,7 +20,6 @@ extern "C" {
 }
 
 #include <cstdint>
-#include <functional>
 #include <unordered_map>
 #include <vector>
 
diff --git a/src/Chemistry/SurrogateModels/InterpolationModule.cpp b/src/Chemistry/SurrogateModels/InterpolationModule.cpp
index 26bfd4ab0..7dc64919b 100644
--- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp
+++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp
@@ -1,8 +1,8 @@
 //  Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
 #include "Interpolation.hpp"
 
-#include "../../DataStructures/DataStructures.hpp"
 #include "DHT_Wrapper.hpp"
+#include "DataStructures/NamedVector.hpp"
 #include "HashFunctions.hpp"
 #include "LookupKey.hpp"
 #include "Rounding.hpp"
diff --git a/src/Chemistry/WorkerFunctions.cpp b/src/Chemistry/WorkerFunctions.cpp
index 8cec96d62..6f4f81c6d 100644
--- a/src/Chemistry/WorkerFunctions.cpp
+++ b/src/Chemistry/WorkerFunctions.cpp
@@ -2,10 +2,8 @@
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include <IrmResult.h>
-#include <algorithm>
-#include <cassert>
-#include <cmath>
+#include "Chemistry/ChemistryDefs.hpp"
+
 #include <cstddef>
 #include <cstdint>
 #include <iomanip>
@@ -18,19 +16,6 @@
 
 namespace poet {
 
-struct WorkPackage {
-  std::size_t size;
-  std::vector<std::vector<double>> input;
-  std::vector<std::vector<double>> output;
-  std::vector<std::uint8_t> mapping;
-
-  WorkPackage(size_t _size) : size(_size) {
-    input.resize(size);
-    output.resize(size);
-    mapping.resize(size, CHEM_PQC);
-  }
-};
-
 inline std::string get_string(int root, MPI_Comm communicator) {
   int count;
   MPI_Bcast(&count, 1, MPI_INT, root, communicator);
@@ -59,13 +44,8 @@ void poet::ChemistryModule::WorkerLoop() {
     PropagateFunctionType(func_type);
 
     switch (func_type) {
-    case CHEM_INIT: {
-      RunInitFile(get_string(0, this->group_comm));
-      break;
-    }
-    case CHEM_INIT_SPECIES: {
-      Field dummy;
-      initializeField(dummy);
+    case CHEM_FIELD_INIT: {
+      ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
       break;
     }
     case CHEM_WORK_LOOP: {
@@ -191,9 +171,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
 
   phreeqc_time_start = MPI_Wtime();
 
-  if (WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt) != IRM_OK) {
-    std::cerr << "Phreeqc error" << std::endl;
-  };
+  WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
 
   phreeqc_time_end = MPI_Wtime();
 
@@ -299,42 +277,65 @@ void poet::ChemistryModule::WorkerReadDHTDump(
   }
 }
 
-IRM_RESULT
-poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
-                                            double dSimTime, double dTimestep) {
+void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
+                                                 double dSimTime,
+                                                 double dTimestep) {
   // check if we actually need to start phreeqc
-  std::vector<std::uint32_t> pqc_mapping;
+  bool queued_cell = false;
 
-  for (std::size_t i = 0; i < work_package.size; i++) {
-    if (work_package.mapping[i] == CHEM_PQC) {
-      pqc_mapping.push_back(i);
+  for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
+    if (work_package.mapping[wp_id] != CHEM_PQC) {
+      continue;
+    }
+
+    const auto &input = work_package.input[wp_id];
+    const auto &sol_id = input[0];
+
+    auto &phreeqc_instance = this->phreeqc_instances[sol_id];
+    work_package.output[wp_id] = work_package.input[wp_id];
+
+    work_package.input[wp_id].erase(work_package.input[wp_id].begin());
+
+    // remove NaNs from the input
+    work_package.input[wp_id].erase(
+        std::remove_if(work_package.input[wp_id].begin(),
+                       work_package.input[wp_id].end(),
+                       [](double d) { return std::isnan(d); }),
+        work_package.input[wp_id].end());
+
+    phreeqc_instance->queueCell(work_package.input[wp_id]);
+    queued_cell = true;
+  }
+
+  if (!queued_cell) {
+    return;
+  }
+
+  std::map<int, std::vector<std::size_t>> zone_mapping;
+
+  // 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++;
     }
   }
-
-  if (pqc_mapping.empty()) {
-    return IRM_OK;
-  }
-
-  IRM_RESULT result;
-  this->PhreeqcRM::setPOETMapping(pqc_mapping);
-  this->setDumpedField(work_package.input);
-
-  this->PhreeqcRM::SetTime(dSimTime);
-  this->PhreeqcRM::SetTimeStep(dTimestep);
-
-  result = this->PhreeqcRM::RunCells();
-
-  std::vector<std::vector<double>> output_tmp(work_package.size);
-
-  this->getDumpedField(output_tmp);
-
-  for (std::size_t i = 0; i < work_package.size; i++) {
-    if (work_package.mapping[i] == CHEM_PQC) {
-      work_package.output[i] = output_tmp[i];
-    }
-  }
-
-  return result;
 }
 
 void poet::ChemistryModule::WorkerPerfToMaster(int type,
diff --git a/src/Chemistry/enums.hpp b/src/Chemistry/enums.hpp
deleted file mode 100644
index 11c3c1281..000000000
--- a/src/Chemistry/enums.hpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#ifndef ENUMS_H_
-#define ENUMS_H_
-
-namespace poet {
-enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
-
-enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
-} // namespace poet
-
-#endif // ENUMS_H_
diff --git a/src/DataStructures/Field.hpp b/src/DataStructures/Field.hpp
index 2f2269653..90119f2a6 100644
--- a/src/DataStructures/Field.hpp
+++ b/src/DataStructures/Field.hpp
@@ -4,6 +4,7 @@
 #include <Rcpp.h>
 
 #include <cstddef>
+#include <cstdint>
 #include <string>
 #include <vector>
 
diff --git a/src/Init/ChemistryInit.cpp b/src/Init/ChemistryInit.cpp
index 2975eec2c..9b86d5c41 100644
--- a/src/Init/ChemistryInit.cpp
+++ b/src/Init/ChemistryInit.cpp
@@ -4,8 +4,6 @@ namespace poet {
 InitialList::ChemistryInit InitialList::getChemistryInit() const {
   ChemistryInit chem_init;
 
-  chem_init.initial_grid = Field(initial_grid);
-
   chem_init.total_grid_cells = this->n_cols * this->n_rows;
 
   chem_init.database = database;
diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp
index 2b9d030ab..2ec7fab33 100644
--- a/src/Init/GridInit.cpp
+++ b/src/Init/GridInit.cpp
@@ -39,9 +39,9 @@ static Rcpp::NumericMatrix pqcScriptToGrid(IPhreeqcPOET &phreeqc, RInside &R) {
 
 static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat,
                                    const Rcpp::NumericMatrix &grid) {
-  Rcpp::Function pqc_to_grid("pqc_to_grid");
+  Rcpp::Function pqc_to_grid_R("pqc_to_grid");
 
-  return pqc_to_grid(mat, grid);
+  return pqc_to_grid_R(mat, grid);
 }
 
 static inline std::map<int, std::string>
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index 86c4ded77..cb11bd3d6 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -160,8 +160,6 @@ private:
 
 public:
   struct ChemistryInit {
-    Field initial_grid;
-
     uint32_t total_grid_cells;
 
     std::string database;
diff --git a/src/initializer.cpp b/src/initializer.cpp
index 1f95babed..eacc45a72 100644
--- a/src/initializer.cpp
+++ b/src/initializer.cpp
@@ -18,7 +18,9 @@ int main(int argc, char **argv) {
   R.parseEvalQ(init_r_library);
 
   const std::string script = argv[1];
-  R.parseEvalQ("source('" + script + "')");
+  Rcpp::Function source_R("source");
+
+  source_R(script);
 
   Rcpp::List setup = R["setup"];
 
diff --git a/src/poet.cpp b/src/poet.cpp
index fa7819460..f162b7820 100644
--- a/src/poet.cpp
+++ b/src/poet.cpp
@@ -23,6 +23,7 @@
 #include "Base/Macros.hpp"
 #include "Base/RInsidePOET.hpp"
 // #include "Chemistry/ChemistryModule.hpp"
+#include "Chemistry/ChemistryModule.hpp"
 #include "DataStructures/Field.hpp"
 #include "Init/InitialList.hpp"
 #include "Transport/DiffusionModule.hpp"
@@ -48,41 +49,16 @@ using namespace Rcpp;
 
 static int MY_RANK = 0;
 
-// poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
-//   std::unordered_map<std::string, double> out_map;
-//   vector<string> col_names = Rcpp::as<vector<string>>(df.names());
-
-//   for (const auto &name : col_names) {
-//     double val = df[name.c_str()];
-//     out_map.insert({name, val});
-//   }
-
-//   return out_map;
-// }
-
 // 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) {
+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");
+
+  R["TMP"] = chem.asSEXP();
   R.parseEval("mysetup$state_C <- TMP");
-
-  // R["TMP"] = Rcpp::wrap(trans.AsVector());
-  // R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
-  // R.parseEval(std::string(
-  //     "mysetup$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(
-  //     "mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
-  //     std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
 }
 
 enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
@@ -102,11 +78,12 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
     return ParseRet::PARSER_HELP;
   }
   // if positional arguments are missing
-  else if (!cmdl(2)) {
+  if (!cmdl(3)) {
     if (MY_RANK == 0) {
-      ERRMSG("POET needs 2 positional arguments: ");
-      ERRMSG("1) the R script defining your simulation and");
-      ERRMSG("2) the directory prefix where to save results and profiling");
+      ERRMSG("POET needs 3 positional arguments: ");
+      ERRMSG("1) the R script defining your simulation runtime.");
+      ERRMSG("2) the initial .rds file generated by poet_init.");
+      ERRMSG("3) the directory prefix where to save results and profiling");
     }
     return ParseRet::PARSER_ERROR;
   }
@@ -203,8 +180,8 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
   std::string runtime_file;
   std::string out_dir;
 
-  cmdl(1) >> init_file;
-  cmdl(2) >> runtime_file;
+  cmdl(1) >> runtime_file;
+  cmdl(2) >> init_file;
   cmdl(3) >> out_dir;
 
   // chem_params.dht_outdir = out_dir;
@@ -300,7 +277,7 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
 // }
 
 static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
-                            DiffusionModule &diffusion) {
+                            DiffusionModule &diffusion, ChemistryModule &chem) {
 
   /* Iteration Count is dynamic, retrieving value from R (is only needed by
    * master for the following loop) */
@@ -320,7 +297,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
   // chem.initializeField(diffusion.getField());
 
   // if (params.getNumParams().print_progressbar) {
-  //   chem.setProgressBarPrintout(true);
+  chem.setProgressBarPrintout(true);
   // }
 
   /* SIMULATION LOOP */
@@ -332,7 +309,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
     // cout << "CPP: Evaluating next time step" << endl;
     // R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
 
-    double dt = params.timesteps[iter - 1];
+    const double &dt = params.timesteps[iter - 1];
 
     //  cout << "CPP: Next time step is " << dt << "[s]" << endl;
     MSG("Next time step is " + std::to_string(dt) + " [s]");
@@ -347,14 +324,15 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
     // TODO: transport to diffusion
     diffusion.simulate(dt);
 
-    // chem.getField().update(diffusion.getField());
+    chem.GetField().update(diffusion.getField());
+
+    // chem.getfield().update(diffusion.getfield());
 
     MSG("Chemistry step");
 
-    // chem.SetTimeStep(dt);
-    // chem.RunCells();
+    chem.simulate(dt);
 
-    writeFieldsToR(R, diffusion.getField());
+    writeFieldsToR(R, diffusion.getField(), chem.GetField());
     // diffusion.getField().update(chem.GetField());
 
     R["req_dt"] = dt;
@@ -431,7 +409,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
   //   R.parseEvalQ("profiling$interp_cached <- interp_cached");
   // }
 
-  // chem.MasterLoopBreak();
+  chem.MasterLoopBreak();
 
   return dSimTime;
 }
@@ -453,37 +431,6 @@ int main(int argc, char *argv[]) {
     MSG("Running POET version " + std::string(poet_version));
   }
 
-  if (world_rank > 0) {
-    {
-      // SimParams params(world_rank, world_size);
-      // int pret = params.parseFromCmdl(argv, R);
-
-      // if (pret == poet::PARSER_ERROR) {
-      //   MPI_Finalize();
-      //   return EXIT_FAILURE;
-      // } else if (pret == poet::PARSER_HELP) {
-      //   MPI_Finalize();
-      //   return EXIT_SUCCESS;
-      // }
-
-      // // ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
-      // ChemistryModule worker = poet::ChemistryModule::createWorker(
-      //     MPI_COMM_WORLD, params.getChemParams());
-      // set_chem_parameters(worker, worker.GetWPSize(),
-      //                     params.getChemParams().database_path);
-
-      // worker.WorkerLoop();
-    }
-
-    MPI_Barrier(MPI_COMM_WORLD);
-
-    MSG("finished, cleanup of process " + std::to_string(world_rank));
-
-    MPI_Finalize();
-
-    return EXIT_SUCCESS;
-  }
-
   /*Loading Dependencies*/
   // TODO: kann raus
   R.parseEvalQ(kin_r_library);
@@ -499,8 +446,25 @@ int main(int argc, char *argv[]) {
     break;
   }
 
+  InitialList init_list(R);
+  init_list.importList(run_params.init_params);
+
   MSG("RInside initialized on process " + std::to_string(world_rank));
 
+  if (world_rank > 0) {
+
+    ChemistryModule worker(1, init_list.getChemistryInit(), MPI_COMM_WORLD);
+    worker.WorkerLoop();
+
+    MPI_Barrier(MPI_COMM_WORLD);
+
+    MSG("finished, cleanup of process " + std::to_string(world_rank));
+
+    MPI_Finalize();
+
+    return EXIT_SUCCESS;
+  }
+
   // R.parseEvalQ("mysetup <- setup");
   // // if (world_rank == 0) { // get timestep vector from
   // // grid_init function ... //
@@ -521,13 +485,14 @@ int main(int argc, char *argv[]) {
 
   // MPI_Barrier(MPI_COMM_WORLD);
 
-  InitialList init_list(R);
-  init_list.importList(run_params.init_params);
-
   DiffusionModule diffusion(init_list.getDiffusionInit(),
                             init_list.getInitialGrid());
 
-  dSimTime = RunMasterLoop(R, run_params, diffusion);
+  ChemistryModule chemistry(1, init_list.getChemistryInit(), MPI_COMM_WORLD);
+
+  chemistry.masterSetField(init_list.getInitialGrid());
+
+  dSimTime = RunMasterLoop(R, run_params, diffusion, chemistry);
 
   MSG("finished simulation loop");