diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R
index 1e7f803ed..d4d04f9e2 100644
--- a/R_lib/kin_r_library.R
+++ b/R_lib/kin_r_library.R
@@ -28,7 +28,7 @@ master_init <- function(setup) {
 
     ## Setup the directory where we will store the results
     verb <- FALSE
-    if (local_rank == 0) {
+    # if (local_rank == 0) {
         verb <- TRUE ## verbosity loading MUFITS results
         if (!dir.exists(fileout)) {
             dir.create(fileout)
@@ -41,25 +41,26 @@ master_init <- function(setup) {
         } else {
             msgm("store_result is ", store_result)
         }
-    } else {
-        
-    }
+    # } else {
+
+    # }
 
     setup$iter <- 1
-    setup$maxiter <- setup$iterations
     setup$timesteps <- setup$timesteps
+    setup$maxiter <- length(setup$timesteps)
+    setup$iterations <- setup$maxiter
     setup$simulation_time <- 0
 
     if (is.null(setup[["store_result"]])) {
         setup$store_result <- TRUE
     }
-    
+
     if (setup$store_result) {
         if (is.null(setup[["out_save"]])) {
             setup$out_save <- seq(1, setup$iterations)
         }
     }
-    
+
     return(setup)
 }
 
@@ -67,31 +68,33 @@ master_init <- function(setup) {
 ## calculated time step if store_result is TRUE and increments the
 ## iteration counter
 master_iteration_end <- function(setup) {
-    iter <- setup$iter
+    # iter <- setup$iter
+    # print(iter)
     ## max digits for iterations
-    dgts <-  as.integer(ceiling(log10(setup$iterations + 1)))
-    ## string format to use in sprintf 
+    iter <- 1
+    dgts <- as.integer(ceiling(log10(1)))
+    ## string format to use in sprintf
     fmt <- paste0("%0", dgts, "d")
-    
+
     ## Write on disk state_T and state_C after every iteration
     ## comprised in setup$out_save
-    if (setup$store_result) {
-        if (iter %in% setup$out_save) {
-            nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds")
-            info <- list(
-                tr_req_dt = as.integer(setup$req_dt)
-                ## tr_allow_dt = setup$allowed_dt,
-                ## tr_inniter = as.integer(setup$inniter)
-            )
-            saveRDS(list(
-                T = setup$state_T, C = setup$state_C,
-                simtime = as.integer(setup$simtime),
-                tr_info = info
-            ), file = nameout)
-            msgm("results stored in <", nameout, ">")
-        }
-    }
-    msgm("done iteration", iter, "/", setup$maxiter)
+    # if (setup$store_result) {
+    #     if (iter %in% setup$out_save) {
+    nameout <- paste0(fileout, "/iter_", sprintf(fmt = fmt, iter), ".rds")
+    info <- list(
+        tr_req_dt = as.integer(1)
+        ## tr_allow_dt = setup$allowed_dt,
+        ## tr_inniter = as.integer(setup$inniter)
+    )
+    saveRDS(list(
+        T = setup$state_T, C = setup$state_C,
+        simtime = as.integer(0),
+        tr_info = info
+    ), file = nameout)
+    msgm("results stored in <", nameout, ">")
+    #     }
+    # }
+    msgm("done iteration", iter, "/", 1)
     setup$iter <- setup$iter + 1
     return(setup)
 }
@@ -104,41 +107,41 @@ slave_chemistry <- function(setup, data) {
     immobile <- setup$immobile
     kin <- setup$kin
     ann <- setup$ann
-    
+
     iter <- setup$iter
     timesteps <- setup$timesteps
     dt <- timesteps[iter]
-    
+
     state_T <- data ## not the global field, but the work-package
-    
+
     ## treat special H+/pH, e-/pe cases
     state_T <- RedModRphree::Act2pH(state_T)
-    
+
     ## reduction of the problem
     if (setup$reduce) {
         reduced <- ReduceStateOmit(state_T, omit = setup$ann)
     } else {
         reduced <- state_T
     }
-    
+
     ## form the PHREEQC input script for the current work package
     inplist <- SplitMultiKin(
         data = reduced, procs = 1, base = base, first = first,
         ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt
     )
-    
+
     ## if (local_rank==1 & iter==1)
     ##         RPhreeWriteInp("FirstInp", inplist)
-    
+
     tmpC <- RunPQC(inplist, procs = 1, second = TRUE)
-    
+
     ## recompose after the reduction
     if (setup$reduce) {
         state_C <- RecomposeState(tmpC, reduced)
     } else {
         state_C <- tmpC
     }
-    
+
     ## the next line is needed since we don't need all columns of
     ## PHREEQC output
     return(state_C[, prop])
@@ -147,40 +150,40 @@ slave_chemistry <- function(setup, data) {
 ## This function, called by master
 master_chemistry <- function(setup, data) {
     state_T <- setup$state_T
-    
+
     msgm(" chemistry iteration", setup$iter)
-    
+
     ## treat special H+/pH, e-/pe cases
     state_T <- RedModRphree::Act2pH(state_T)
-    
+
     ## reduction of the problem
     if (setup$reduce) {
         reduced <- ReduceStateOmit(state_T, omit = setup$ann)
     } else {
         reduced <- state_T
     }
-    
+
     ## inject data from workers
     res_C <- data
-    
+
     rownames(res_C) <- NULL
-    
+
     ## print(res_C)
-    
+
     if (nrow(res_C) > nrow(reduced)) {
         res_C <- res_C[seq(2, nrow(res_C), by = 2), ]
     }
-    
+
     ## recompose after the reduction
     if (setup$reduce) {
         state_C <- RecomposeState(res_C, reduced)
     } else {
         state_C <- res_C
     }
-    
+
     setup$state_C <- state_C
     setup$reduced <- reduced
-    
+
     return(setup)
 }
 
@@ -188,7 +191,7 @@ master_chemistry <- function(setup, data) {
 ## Adapted version for "reduction"
 ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
     require(mgcv)
-    
+
     rem <- colnames(data)
     if (is.list(omit)) {
         indomi <- match(names(omit), colnames(data))
@@ -196,13 +199,13 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
     } else {
         datao <- data
     }
-    
+
     datao <- signif(datao, sign)
     red <- mgcv::uniquecombs(datao)
     inds <- attr(red, "index")
     now <- ncol(red)
-    
-    
+
+
     ## reattach the omitted column(s)
     ## FIXME: control if more than one ann is present
     if (is.list(omit)) {
@@ -224,11 +227,11 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
 ## Attach the name of the calling function to the message displayed on
 ## R's stdout
 msgm <- function(...) {
-    if (local_rank == 0) {
-        fname <- as.list(sys.call(-1))[[1]]
-        prefix <- paste0("R: ", fname, " ::")
-        cat(paste(prefix, ..., "\n"))
-    }
+    # if (local_rank == 0) {
+    fname <- as.list(sys.call(-1))[[1]]
+    prefix <- paste0("R: ", fname, " ::")
+    cat(paste(prefix, ..., "\n"))
+    # }
     invisible()
 }
 
@@ -236,18 +239,17 @@ msgm <- function(...) {
 ## Function called by master R process to store on disk all relevant
 ## parameters for the simulation
 StoreSetup <- function(setup) {
-
     to_store <- vector(mode = "list", length = 4)
     ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
     names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
-    
+
     ## read the setup R file, which is sourced in kin.cpp
     tmpbuff <- file(filesim, "r")
     setupfile <- readLines(tmpbuff)
     close.connection(tmpbuff)
-    
+
     to_store$Sim <- setupfile
-    
+
     ## to_store$Flow <- list(
     ##     snapshots  = setup$snapshots,
     ##     gridfile   = setup$gridfile,
@@ -282,22 +284,22 @@ StoreSetup <- function(setup) {
         to_store$DHT <- FALSE
     }
 
-  if (dht_enabled) {
-    to_store$DHT <- list(
-      enabled   = dht_enabled,
-      log       = dht_log
-      #signif    = dht_final_signif,
-      #proptype  = dht_final_proptype
-    )
-  } else {
-    to_store$DHT <- FALSE
-  }
+    if (dht_enabled) {
+        to_store$DHT <- list(
+            enabled   = dht_enabled,
+            log       = dht_log
+            # signif    = dht_final_signif,
+            # proptype  = dht_final_proptype
+        )
+    } else {
+        to_store$DHT <- FALSE
+    }
 
-  saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
-  msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
+    saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
+    msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
 }
 
 GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
-    ids <- rep(1:n_packages, times=package_size, each = 1)[1:len]
+    ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
     return(as.integer(table(ids)))
 }
diff --git a/bench/barite/het/runtime.R b/bench/barite/het/runtime.R
new file mode 100644
index 000000000..4dea41fdf
--- /dev/null
+++ b/bench/barite/het/runtime.R
@@ -0,0 +1,4 @@
+list(
+    timesteps = c(1),
+    store_result = TRUE
+)
\ No newline at end of file
diff --git a/ext/tug b/ext/tug
index 25855da6b..8c0687a6c 160000
--- a/ext/tug
+++ b/ext/tug
@@ -1 +1 @@
-Subproject commit 25855da6b2930559b542bbadb16299932332d6a3
+Subproject commit 8c0687a6cd4a10a79c7a554083a35eda11cc55f0
diff --git a/src/Base/Grid.cpp b/src/Base/Grid.cpp
deleted file mode 100644
index e51b114ba..000000000
--- a/src/Base/Grid.cpp
+++ /dev/null
@@ -1,274 +0,0 @@
-/*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-#include "Grid.hpp"
-#include "SimParams.hpp"
-
-#include <RInside.h>
-#include <Rcpp.h>
-#include <algorithm>
-#include <cassert>
-#include <cstdint>
-#include <new>
-#include <numeric>
-#include <stdexcept>
-#include <string>
-#include <vector>
-
-using namespace poet;
-using namespace Rcpp;
-
-enum { INIT_SCRATCH, INIT_PHREEQC, INIT_RDS };
-
-static inline int8_t get_type_id(std::string type) {
-  if (type == "scratch") {
-    return INIT_SCRATCH;
-  }
-  if (type == "phreeqc") {
-    return INIT_PHREEQC;
-  }
-  if (type == "rds") {
-    return INIT_RDS;
-  }
-  return -1;
-}
-
-void poet::Grid::InitModuleFromParams(const poet::GridParams &grid_args) {
-  assert(grid_args.n_cells.size() == grid_args.s_cells.size());
-
-  this->SetGridDimension(grid_args.n_cells.size());
-
-  std::copy(grid_args.n_cells.begin(), grid_args.n_cells.end(),
-            this->n_cells.begin());
-  std::copy(grid_args.s_cells.begin(), grid_args.s_cells.end(),
-            this->grid_size.begin());
-
-  auto init_type = get_type_id(grid_args.type);
-  assert(init_type >= 0);
-
-  switch (init_type) {
-  case INIT_SCRATCH:
-    this->SetPropNames(grid_args.props);
-    this->InitGridFromScratch(grid_args.init_df);
-    break;
-  case INIT_RDS:
-    this->SetPropNames(grid_args.props);
-    this->InitGridFromRDS();
-    break;
-  case INIT_PHREEQC:
-    this->InitGridFromPhreeqc();
-  }
-}
-
-void poet::Grid::SetGridDimension(uint8_t dim) {
-  assert(dim < MAX_DIM + 1);
-
-  this->dim = dim;
-}
-
-void poet::Grid::SetGridCellCount(uint32_t n_x, uint32_t n_y, uint32_t n_z) {
-  assert(this->dim < MAX_DIM + 1);
-
-  this->n_cells[0] = n_x;
-  this->n_cells[1] = n_y;
-  // this->n_cells[2] = n_z;
-}
-
-void poet::Grid::SetGridSize(double s_x, double s_y, double s_z) {
-  assert(this->dim < MAX_DIM + 1);
-
-  this->grid_size[0] = s_x;
-  this->grid_size[1] = s_y;
-  // this->grid_size[2] = s_z;
-}
-
-void poet::Grid::SetPropNames(const std::vector<std::string> &prop_names) {
-  this->prop_names = prop_names;
-}
-
-void poet::Grid::PushbackModuleFlow(const std::string &input,
-                                    const std::string &output) {
-  FlowInputOutputInfo element = {input, output};
-  this->flow_vec.push_back(element);
-}
-
-void poet::Grid::PreModuleFieldCopy(uint32_t tick) {
-  FlowInputOutputInfo curr_element = this->flow_vec.at(tick);
-
-  const std::string input_module_name = curr_element.input_field;
-  StateMemory *out_state = this->GetStatePointer(curr_element.output_field);
-
-  std::vector<std::string> &mod_props = out_state->props;
-  std::vector<double> &mod_field = out_state->mem;
-
-  // copy output of another module  as input for another module
-  for (uint32_t i = 0; i < mod_props.size(); i++) {
-    try {
-      std::vector<double> t_prop_vec =
-          this->GetSpeciesByName(mod_props[i], input_module_name);
-
-      std::copy(t_prop_vec.begin(), t_prop_vec.end(),
-                mod_field.begin() + (i * this->GetTotalCellCount()));
-    } catch (...) {
-      // TODO: there might be something wrong ...
-      continue;
-    }
-  }
-}
-
-Grid::Grid() {
-  this->n_cells.fill(0);
-  this->grid_size.fill(0);
-};
-
-Grid::~Grid() {
-  for (auto &[key, val] : this->state_register) {
-    delete val;
-  }
-}
-
-void poet::Grid::InitGridFromScratch(const Rcpp::DataFrame &init_cell) {
-  const uint32_t vec_size = this->GetTotalCellCount();
-
-  StateMemory *curr_state =
-      this->RegisterState(poet::GRID_MODULE_NAME, this->prop_names);
-  curr_state->props = this->prop_names;
-
-  std::vector<double> &curr_field = curr_state->mem;
-
-  for (auto const &prop : this->prop_names) {
-    // size of the vector shall be 1
-    double prop_val = ((Rcpp::NumericVector)init_cell[prop.c_str()])[0];
-    std::vector<double> prop_vec(vec_size, prop_val);
-
-    curr_field.insert(curr_field.end(), prop_vec.begin(), prop_vec.end());
-  }
-
-  this->grid_init = curr_state;
-}
-void poet::Grid::InitGridFromRDS() {
-  // TODO
-}
-void poet::Grid::InitGridFromPhreeqc() {
-  // TODO
-}
-
-poet::StateMemory *Grid::RegisterState(std::string module_name,
-                                       std::vector<std::string> props) {
-  poet::StateMemory *mem = new poet::StateMemory;
-
-  mem->props = props;
-
-  const auto [it, success] = this->state_register.insert({module_name, mem});
-
-  if (!success) {
-    delete mem;
-    throw std::bad_alloc();
-  }
-
-  return mem;
-}
-
-poet::StateMemory *Grid::GetStatePointer(std::string module_name) {
-  const auto it = this->state_register.find(module_name);
-
-  if (it == this->state_register.end()) {
-    throw std::range_error(
-        std::string("Requested module " + module_name + " was not found"));
-  }
-
-  return it->second;
-}
-
-auto Grid::GetGridSize(uint8_t direction) const -> uint32_t {
-  return this->grid_size.at(direction);
-}
-
-auto Grid::GetGridDimension() const -> uint8_t { return this->dim; }
-
-auto Grid::GetTotalCellCount() const -> uint32_t {
-  uint32_t sum = 1;
-  for (auto const &dim : this->n_cells) {
-    sum *= (dim != 0 ? dim : 1);
-  }
-
-  return sum;
-}
-
-auto Grid::GetGridCellsCount(uint8_t direction) const -> uint32_t {
-  return this->n_cells.at(direction);
-}
-
-auto poet::Grid::GetInitialGrid() const -> StateMemory * {
-  return this->grid_init;
-}
-
-auto Grid::GetSpeciesCount() const -> uint32_t {
-  return this->prop_names.size();
-}
-
-auto Grid::GetPropNames() const -> std::vector<std::string> {
-  return this->prop_names;
-}
-
-auto poet::Grid::GetSpeciesByName(std::string name,
-                                  std::string module_name) const
-    -> std::vector<double> {
-
-  auto const it = this->state_register.find(module_name);
-
-  if (it == this->state_register.end()) {
-    throw std::range_error(
-        std::string("Requested module " + module_name + " was not found"));
-  }
-
-  poet::StateMemory const *module_memory = it->second;
-
-  auto const &prop_vector = module_memory->props;
-  auto const find_res = std::find(prop_vector.begin(), prop_vector.end(), name);
-  if (find_res == prop_vector.end()) {
-    throw std::range_error(std::string(
-        "Species " + name + " was not found for module " + module_name));
-  }
-  uint32_t prop_index = find_res - prop_vector.begin();
-
-  uint32_t begin_vec = prop_index * this->GetTotalCellCount(),
-           end_vec = ((prop_index + 1) * this->GetTotalCellCount());
-
-  return std::vector<double>(module_memory->mem.begin() + begin_vec,
-                             module_memory->mem.begin() + end_vec);
-}
-
-// HACK: Helper function
-void poet::Grid::WriteFieldsToR(RInside &R) {
-
-  for (auto const &elem : this->state_register) {
-    std::string field_name = elem.first;
-    StateMemory *field = elem.second;
-
-    const uint32_t n_cells_per_prop = field->mem.size() / field->props.size();
-
-    R["TMP"] = Rcpp::wrap(field->mem);
-    R["TMP_PROPS"] = Rcpp::wrap(field->props);
-    R.parseEval(std::string(
-        "mysetup$" + field_name + " <- setNames(data.frame(matrix(TMP, nrow=" +
-        std::to_string(n_cells_per_prop) + ")), TMP_PROPS)"));
-  }
-}
diff --git a/src/Base/Grid.hpp b/src/Base/Grid.hpp
deleted file mode 100644
index deec5c04a..000000000
--- a/src/Base/Grid.hpp
+++ /dev/null
@@ -1,117 +0,0 @@
-/*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-#ifndef GRID_H
-#define GRID_H
-
-#include "SimParams.hpp"
-
-#include <RInside.h>
-#include <Rcpp.h>
-#include <array>
-#include <cstddef>
-#include <cstdint>
-#include <list>
-#include <map>
-#include <string>
-#include <vector>
-
-#define MAX_DIM 2
-
-namespace poet {
-
-enum { GRID_X_DIR, GRID_Y_DIR, GRID_Z_DIR };
-
-using StateMemory = struct s_StateMemory {
-  std::vector<double> mem;
-  std::vector<std::string> props;
-};
-
-using FlowInputOutputInfo = struct inOut_info {
-  std::string input_field;
-  std::string output_field;
-};
-
-constexpr const char *GRID_MODULE_NAME = "grid_init";
-
-/**
- * @brief Class describing the grid
- *
- * Providing methods to shuffle and unshuffle grid (for the master) as also to
- * import and export a work package (for worker).
- *
- * @todo find better abstraction
- *
- */
-class Grid {
-
-public:
-  Grid();
-
-  ~Grid();
-
-  void InitModuleFromParams(const poet::GridParams &grid_args);
-
-  void SetGridDimension(uint8_t dim);
-  void SetGridCellCount(uint32_t n_x, uint32_t n_y = 0, uint32_t n_z = 0);
-  void SetGridSize(double s_x, double s_y = 0., double s_z = 0.);
-  void SetPropNames(const std::vector<std::string> &prop_names);
-
-  void PushbackModuleFlow(const std::string &input, const std::string &output);
-  void PreModuleFieldCopy(uint32_t tick);
-
-  void InitGridFromScratch(const Rcpp::DataFrame &init_cell);
-  void InitGridFromRDS();
-  void InitGridFromPhreeqc();
-
-  auto GetGridDimension() const -> uint8_t;
-  auto GetTotalCellCount() const -> uint32_t;
-  auto GetGridCellsCount(uint8_t direction) const -> uint32_t;
-  auto GetGridSize(uint8_t direction) const -> uint32_t;
-
-  auto RegisterState(std::string module_name, std::vector<std::string> props)
-      -> StateMemory *;
-  auto GetStatePointer(std::string module_name) -> StateMemory *;
-
-  auto GetInitialGrid() const -> StateMemory *;
-
-  auto GetSpeciesCount() const -> uint32_t;
-  auto GetPropNames() const -> std::vector<std::string>;
-
-  auto GetSpeciesByName(std::string name,
-                        std::string module_name = poet::GRID_MODULE_NAME) const
-      -> std::vector<double>;
-
-  void WriteFieldsToR(RInside &R);
-
-private:
-  std::vector<FlowInputOutputInfo> flow_vec;
-
-  std::uint8_t dim = 0;
-  std::array<double, MAX_DIM> grid_size;
-  std::array<std::uint32_t, MAX_DIM> n_cells;
-
-  std::map<std::string, StateMemory *> state_register;
-  StateMemory *grid_init = std::nullptr_t();
-
-  std::vector<std::string> prop_names;
-};
-} // namespace poet
-#endif // GRID_H
diff --git a/src/Base/SimParams.cpp b/src/Base/SimParams.cpp
index f711e796a..3ae76dd65 100644
--- a/src/Base/SimParams.cpp
+++ b/src/Base/SimParams.cpp
@@ -1,10 +1,11 @@
-//  Time-stamp: "Last modified 2023-08-15 12:12:36 delucia"
-
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
 ** Potsdam)
 **
-** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
+** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
+**
+** Copyright (C) 2018-2022 Marco De Lucia (GFZ Potsdam), Max Luebke (University
+** of Potsdam)
 **
 ** POET is free software; you can redistribute it and/or modify it under the
 ** terms of the GNU General Public License as published by the Free Software
@@ -22,104 +23,33 @@
 
 #include "SimParams.hpp"
 
-#include "../Chemistry/enums.hpp"
+#include "Base/Macros.hpp"
 
-#include <algorithm>
-#include <cassert>
 #include <cstdint>
 
-#include <RInside.h>
-#include <Rcpp.h>
-#include <algorithm>
-#include <cassert>
 #include <cstdint>
-#include <iostream>
 #include <string>
-#include <string_view>
 #include <vector>
 
+#include "Base/RInsidePOET.hpp"
+#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
+
 using namespace poet;
-using namespace std;
-using namespace Rcpp;
 
-poet::GridParams::s_GridParams(RInside &R) {
-  auto tmp_n_cells =
-      Rcpp::as<std::vector<uint32_t>>(R.parseEval("mysetup$grid$n_cells"));
-  assert(tmp_n_cells.size() < 3);
+RuntimeParameters::RuntimeParameters(RInsidePOET &R, char *argv[],
+                                     int world_rank_)
+    : world_rank(world_rank_) {
 
-  this->dim = tmp_n_cells.size();
-
-  std::copy(tmp_n_cells.begin(), tmp_n_cells.end(), this->n_cells.begin());
-
-  auto tmp_s_cells =
-      Rcpp::as<std::vector<double>>(R.parseEval("mysetup$grid$s_cells"));
-
-  assert(tmp_s_cells.size() == this->dim);
-
-  std::copy(tmp_s_cells.begin(), tmp_s_cells.end(), this->s_cells.begin());
-
-  this->total_n =
-      (dim == 1 ? this->n_cells[0] : this->n_cells[0] * this->n_cells[1]);
-
-  this->type = Rcpp::as<std::string>(R.parseEval("mysetup$grid$type"));
-}
-
-poet::DiffusionParams::s_DiffusionParams(RInside &R) {
-  this->initial_t =
-      Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$init"));
-  this->alpha =
-      Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
-  if (Rcpp::as<bool>(
-          R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) {
-    this->vecinj_inner =
-        Rcpp::as<Rcpp::List>(R.parseEval("mysetup$diffusion$vecinj_inner"));
-  }
-  this->vecinj =
-      Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
-  this->vecinj_index =
-      Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
-}
-
-void poet::ChemistryParams::initFromR(RInsidePOET &R) {
-  this->database_path =
-      Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
-  this->input_script =
-      Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
-
-  if (R.checkIfExists("dht_species", "mysetup$chemistry")) {
-    this->dht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
-        R.parseEval(("mysetup$chemistry$dht_species")));
-  }
-
-  if (R.checkIfExists("pht_species", "mysetup$chemistry")) {
-    this->pht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
-        R.parseEval(("mysetup$chemistry$pht_species")));
-  }
-  this->hooks.dht_fill =
-      RHookFunction<bool>(R, "mysetup$chemistry$hooks$dht_fill");
-  this->hooks.dht_fuzz =
-      RHookFunction<std::vector<double>>(R, "mysetup$chemistry$hooks$dht_fuzz");
-  this->hooks.interp_pre = RHookFunction<std::vector<std::size_t>>(
-      R, "mysetup$chemistry$hooks$interp_pre_func");
-  this->hooks.interp_post =
-      RHookFunction<bool>(R, "mysetup$chemistry$hooks$interp_post_func");
-}
-
-SimParams::SimParams(int world_rank_, int world_size_) {
-  this->simparams.world_rank = world_rank_;
-  this->simparams.world_size = world_size_;
-}
-
-int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) {
   // initialize argh object
   argh::parser cmdl(argv);
 
   // if user asked for help
   if (cmdl[{"help", "h"}]) {
-    if (simparams.world_rank == 0) {
+    if (this->world_rank == 0) {
       MSG("Todo");
       MSG("See README.md for further information.");
     }
+
     return poet::PARSER_HELP;
   }
   // if positional arguments are missing
@@ -242,11 +172,74 @@ int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) {
   this->chem_params.initFromR(R);
 
   return poet::PARSER_OK;
+};
+
+// poet::GridParams::s_GridParams(RInside &R) {
+//   auto tmp_n_cells =
+//       Rcpp::as<std::vector<uint32_t>>(R.parseEval("mysetup$grid$n_cells"));
+//   assert(tmp_n_cells.size() < 3);
+
+//   this->dim = tmp_n_cells.size();
+
+//   std::copy(tmp_n_cells.begin(), tmp_n_cells.end(), this->n_cells.begin());
+
+//   auto tmp_s_cells =
+//       Rcpp::as<std::vector<double>>(R.parseEval("mysetup$grid$s_cells"));
+
+//   assert(tmp_s_cells.size() == this->dim);
+
+//   std::copy(tmp_s_cells.begin(), tmp_s_cells.end(), this->s_cells.begin());
+
+//   this->total_n =
+//       (dim == 1 ? this->n_cells[0] : this->n_cells[0] * this->n_cells[1]);
+
+//   this->type = Rcpp::as<std::string>(R.parseEval("mysetup$grid$type"));
+// }
+
+// poet::DiffusionParams::s_DiffusionParams(RInside &R) {
+//   this->initial_t =
+//       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$init"));
+//   this->alpha =
+//       Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
+//   if (Rcpp::as<bool>(
+//           R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) {
+//     this->vecinj_inner =
+//         Rcpp::as<Rcpp::List>(R.parseEval("mysetup$diffusion$vecinj_inner"));
+//   }
+//   this->vecinj =
+//       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
+//   this->vecinj_index =
+//       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
+// }
+
+void poet::ChemistryParams::initFromR(RInsidePOET &R) {
+  // this->database_path =
+  //     Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
+  // this->input_script =
+  //     Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
+
+  if (R.checkIfExists("dht_species", "mysetup$chemistry")) {
+    this->dht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
+        R.parseEval(("mysetup$chemistry$dht_species")));
+  }
+
+  if (R.checkIfExists("pht_species", "mysetup$chemistry")) {
+    this->pht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
+        R.parseEval(("mysetup$chemistry$pht_species")));
+  }
+  this->hooks.dht_fill =
+      RHookFunction<bool>(R, "mysetup$chemistry$hooks$dht_fill");
+  this->hooks.dht_fuzz =
+      RHookFunction<std::vector<double>>(R, "mysetup$chemistry$hooks$dht_fuzz");
+  this->hooks.interp_pre = RHookFunction<std::vector<std::size_t>>(
+      R, "mysetup$chemistry$hooks$interp_pre_func");
+  this->hooks.interp_post =
+      RHookFunction<bool>(R, "mysetup$chemistry$hooks$interp_post_func");
 }
 
-void SimParams::initVectorParams(RInside &R) {}
+void RuntimeParameters::initVectorParams(RInside &R) {}
 
-std::list<std::string> SimParams::validateOptions(argh::parser cmdl) {
+std::list<std::string> RuntimeParameters::validateOptions(argh::parser cmdl) {
   /* store all unknown parameters here */
   std::list<std::string> retList;
 
diff --git a/src/Base/SimParams.hpp b/src/Base/SimParams.hpp
index 523b5b105..b4e20f262 100644
--- a/src/Base/SimParams.hpp
+++ b/src/Base/SimParams.hpp
@@ -4,6 +4,9 @@
 **
 ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
 **
+** Copyright (C) 2018-2022 Marco De Lucia (GFZ Potsdam), Max Luebke (University
+** of Potsdam)
+**
 ** POET is free software; you can redistribute it and/or modify it under the
 ** terms of the GNU General Public License as published by the Free Software
 ** Foundation; either version 2 of the License, or (at your option) any later
@@ -18,20 +21,13 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
-#ifndef PARSER_H
-#define PARSER_H
+#pragma once
 
-#include "../DataStructures/DataStructures.hpp"
-#include "Macros.hpp"
+#include "DataStructures/NamedVector.hpp"
 #include "RInsidePOET.hpp"
 
-#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
-
 #include <cstdint>
-#include <optional>
 #include <string>
-#include <string_view>
-#include <unordered_map>
 #include <vector>
 
 #include <RInside.h>
@@ -51,7 +47,8 @@ enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
  * @brief Defining all simulation parameters
  *
  */
-typedef struct {
+struct RuntimeParameters {
+
   /** Count of processes in MPI_COMM_WORLD */
   int world_size;
   /** rank of proces in MPI_COMM_WORLD */
@@ -79,40 +76,40 @@ typedef struct {
   bool print_progressbar;
 
   bool interp_enabled;
-} t_simparams;
-
-using GridParams = struct s_GridParams {
-  std::array<uint32_t, 2> n_cells;
-  std::array<double, 2> s_cells;
-  std::uint8_t dim;
-  std::uint32_t total_n;
-
-  std::string type;
-
-  Rcpp::DataFrame init_df;
-  std::string input_script;
-  std::string database_path;
-
-  std::vector<std::string> props;
-
-  s_GridParams(RInside &R);
 };
 
-using DiffusionParams = struct s_DiffusionParams {
-  Rcpp::DataFrame initial_t;
+// using GridParams = struct s_GridParams {
+//   std::array<uint32_t, 2> n_cells;
+//   std::array<double, 2> s_cells;
+//   std::uint8_t dim;
+//   std::uint32_t total_n;
 
-  Rcpp::NumericVector alpha;
-  Rcpp::List vecinj_inner;
+//   std::string type;
 
-  Rcpp::DataFrame vecinj;
-  Rcpp::DataFrame vecinj_index;
+//   Rcpp::DataFrame init_df;
+//   std::string input_script;
+//   std::string database_path;
 
-  s_DiffusionParams(RInside &R);
-};
+//   std::vector<std::string> props;
+
+//   s_GridParams(RInside &R);
+// };
+
+// using DiffusionParams = struct s_DiffusionParams {
+//   Rcpp::DataFrame initial_t;
+
+//   Rcpp::NumericVector alpha;
+//   Rcpp::List vecinj_inner;
+
+//   Rcpp::DataFrame vecinj;
+//   Rcpp::DataFrame vecinj_index;
+
+//   s_DiffusionParams(RInside &R);
+// };
 
 struct ChemistryParams {
-  std::string database_path;
-  std::string input_script;
+  // std::string database_path;
+  // std::string input_script;
 
   bool use_dht;
   std::uint64_t dht_size;
@@ -137,6 +134,7 @@ struct ChemistryParams {
   void initFromR(RInsidePOET &R);
 };
 
+} // namespace poet
 /**
  * @brief Reads information from program arguments and R runtime
  *
@@ -147,129 +145,136 @@ struct ChemistryParams {
  * Stores and distribute current simulation parameters at any time.
  *
  */
-class SimParams {
-public:
-  /**
-   * @brief Construct a new SimParams object
-   *
-   * With all given parameters a new instance of this class will be created.
-   *
-   * @param world_rank Rank of process inside MPI_COMM_WORLD
-   * @param world_size Size of communicator MPI_COMM_WORLD
-   */
-  SimParams(int world_rank, int world_size);
+// class SimParams {
+// public:
+//   /**
+//    * @brief Construct a new SimParams object
+//    *
+//    * With all given parameters a new instance of this class will be created.
+//    *
+//    * @param world_rank Rank of process inside MPI_COMM_WORLD
+//    * @param world_size Size of communicator MPI_COMM_WORLD
+//    */
+//   SimParams(int world_rank);
 
-  /**
-   * @brief Parse program arguments
-   *
-   * This is done by the argh.h library.
-   *
-   * First, the function will check if there is a flag 'help' or 'h'. If this is
-   * the case a help message is printed and the function will return with
-   * PARSER_HELP.
-   *
-   * Second, if there are not 2 positional arguments an error will be printed to
-   * stderr and the function returns with PARSER_ERROR.
-   *
-   * Then all given program parameters and flags will be read and checked, if
-   * there are known by validateOptions. A list of all unknown options might be
-   * returned, printed out and the function will return with PARSER_ERROR.
-   * Oterhwise the function continuos.
-   *
-   * Now all program arguments will be stored inside t_simparams struct, printed
-   * out and the function returns with PARSER_OK.
-   *
-   * Also, all parsed agruments are distributed to the R runtime.
-   *
-   * @param argv Argument value of the program
-   * @param R Instantiated R runtime
-   * @return int Returns with 0 if no error occured, otherwise value less than 0
-   * is returned.
-   */
-  int parseFromCmdl(char *argv[], RInsidePOET &R);
+//   /**
+//    * @brief Parse program arguments
+//    *
+//    * This is done by the argh.h library.
+//    *
+//    * First, the function will check if there is a flag 'help' or 'h'. If this
+//    is
+//    * the case a help message is printed and the function will return with
+//    * PARSER_HELP.
+//    *
+//    * Second, if there are not 2 positional arguments an error will be printed
+//    to
+//    * stderr and the function returns with PARSER_ERROR.
+//    *
+//    * Then all given program parameters and flags will be read and checked, if
+//    * there are known by validateOptions. A list of all unknown options might
+//    be
+//    * returned, printed out and the function will return with PARSER_ERROR.
+//    * Oterhwise the function continuos.
+//    *
+//    * Now all program arguments will be stored inside t_simparams struct,
+//    printed
+//    * out and the function returns with PARSER_OK.
+//    *
+//    * Also, all parsed agruments are distributed to the R runtime.
+//    *
+//    * @param argv Argument value of the program
+//    * @param R Instantiated R runtime
+//    * @return int Returns with 0 if no error occured, otherwise value less
+//    than 0
+//    * is returned.
+//    */
+//   int parseFromCmdl(char *argv[], RInsidePOET &R);
 
-  /**
-   * @brief Init std::vector values
-   *
-   * This will initialize dht_signif_vector and dht_prop_type_vector internally
-   * depending on whether vectors are defined by R-Simulation file or not.
-   * 'init_chemistry' must be run beforehand.
-   *
-   * @param R R runtime
-   */
-  void initVectorParams(RInside &R);
+//   /**
+//    * @brief Init std::vector values
+//    *
+//    * This will initialize dht_signif_vector and dht_prop_type_vector
+//    internally
+//    * depending on whether vectors are defined by R-Simulation file or not.
+//    * 'init_chemistry' must be run beforehand.
+//    *
+//    * @param R R runtime
+//    */
+//   void initVectorParams(RInside &R);
 
-  /**
-   * @brief Get the numerical params struct
-   *
-   * Returns a struct which contains all numerical or boolean simulation
-   * parameters.
-   *
-   * @return t_simparams Parameter struct
-   */
-  auto getNumParams() const { return this->simparams; };
+//   /**
+//    * @brief Get the numerical params struct
+//    *
+//    * Returns a struct which contains all numerical or boolean simulation
+//    * parameters.
+//    *
+//    * @return t_simparams Parameter struct
+//    */
+//   auto getNumParams() const { return this->simparams; };
 
-  /**
-   * @brief Get the DHT_Signif_Vector
-   *
-   * Returns a vector indicating which significant values are used for each
-   * variable of a grid cell.
-   *
-   * @return std::vector<int> Vector of integers containing information about
-   * significant digits
-   */
-  auto getDHTSignifVector() const { return this->dht_signif_vector; };
+//   /**
+//    * @brief Get the DHT_Signif_Vector
+//    *
+//    * Returns a vector indicating which significant values are used for each
+//    * variable of a grid cell.
+//    *
+//    * @return std::vector<int> Vector of integers containing information about
+//    * significant digits
+//    */
+//   auto getDHTSignifVector() const { return this->dht_signif_vector; };
 
-  auto getPHTSignifVector() const { return this->pht_signif_vector; };
+//   auto getPHTSignifVector() const { return this->pht_signif_vector; };
 
-  auto getPHTBucketSize() const { return this->pht_bucket_size; };
-  auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed; };
+//   auto getPHTBucketSize() const { return this->pht_bucket_size; };
+//   auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed;
+//   };
 
-  /**
-   * @brief Get the filesim name
-   *
-   * Returns a string containing the absolute path to a R file defining the
-   * simulation.
-   *
-   * @return std::string Absolute path to R file
-   */
-  auto getFilesim() const { return this->filesim; };
+//   /**
+//    * @brief Get the filesim name
+//    *
+//    * Returns a string containing the absolute path to a R file defining the
+//    * simulation.
+//    *
+//    * @return std::string Absolute path to R file
+//    */
+//   auto getFilesim() const { return this->filesim; };
 
-  /**
-   * @brief Get the output directory
-   *
-   * Returns the name of an absolute path where all output files should be
-   * stored.
-   *
-   * @return std::string Absolute path to output path
-   */
-  auto getOutDir() const { return this->out_dir; };
+//   /**
+//    * @brief Get the output directory
+//    *
+//    * Returns the name of an absolute path where all output files should be
+//    * stored.
+//    *
+//    * @return std::string Absolute path to output path
+//    */
+//   auto getOutDir() const { return this->out_dir; };
 
-  const auto &getChemParams() const { return chem_params; }
+//   const auto &getChemParams() const { return chem_params; }
 
-private:
-  std::list<std::string> validateOptions(argh::parser cmdl);
+// private:
+//   std::list<std::string> validateOptions(argh::parser cmdl);
 
-  const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
-                                       "interp"};
-  const std::set<std::string> paramlist{
-      "work-package-size", "dht-strategy",
-      "dht-size",          "dht-snaps",
-      "dht-file",          "interp-size",
-      "interp-min",        "interp-bucket-entries"};
+//   const std::set<std::string> flaglist{"ignore-result", "dht", "P",
+//   "progress",
+//                                        "interp"};
+//   const std::set<std::string> paramlist{
+//       "work-package-size", "dht-strategy",
+//       "dht-size",          "dht-snaps",
+//       "dht-file",          "interp-size",
+//       "interp-min",        "interp-bucket-entries"};
 
-  t_simparams simparams;
+//   t_simparams simparams;
 
-  std::vector<uint32_t> dht_signif_vector;
-  std::vector<uint32_t> pht_signif_vector;
+//   std::vector<uint32_t> dht_signif_vector;
+//   std::vector<uint32_t> pht_signif_vector;
 
-  uint32_t pht_bucket_size;
-  uint32_t pht_min_entries_needed;
+//   uint32_t pht_bucket_size;
+//   uint32_t pht_min_entries_needed;
 
-  std::string filesim;
-  std::string out_dir;
+//   std::string filesim;
+//   std::string out_dir;
 
-  ChemistryParams chem_params;
-};
-} // namespace poet
-#endif // PARSER_H
+//   ChemistryParams chem_params;
+// };
+// } // namespace poet
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f593d6b63..2e4430848 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,14 +1,13 @@
-add_library(POETLib OBJECT)
+add_library(POETLib)
 
-target_sources(
-  POETLib 
+target_sources(POETLib 
   PRIVATE
   Init/InitialList.cpp 
   Init/GridInit.cpp 
   Init/DiffusionInit.cpp
   Init/ChemistryInit.cpp
-  PRIVATE 
   DataStructures/Field.cpp
+  Transport/DiffusionModule.cpp
 )
 
 target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
@@ -16,7 +15,7 @@ target_link_libraries(
   POETLib 
   PUBLIC RRuntime
   PUBLIC IPhreeqcPOET
-  PRIVATE tug
+  PUBLIC tug
 )
 
 # add_library(poetlib
@@ -43,7 +42,7 @@ target_link_libraries(
 #   tug
 # )
 
-# target_compile_definitions(poetlib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
+target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
 
 # mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
 # set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
@@ -65,9 +64,9 @@ file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB)
 
 configure_file(poet.hpp.in poet.hpp @ONLY)
 
-# add_executable(poet poet.cpp)
-# target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime)
-# target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
+add_executable(poet poet.cpp)
+target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime)
+target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
 
 add_executable(poet_init initializer.cpp)
 target_link_libraries(poet_init PRIVATE POETLib RRuntime)
diff --git a/src/Init/CMakeLists.txt b/src/Init/CMakeLists.txt
deleted file mode 100644
index 47ce9eb3e..000000000
--- a/src/Init/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-target_sources(POETLib PRIVATE 
-    InitialList.cpp 
-    GridInit.cpp
-    DiffusionInit.cpp
-    ChemistryInit.cpp
-)
\ No newline at end of file
diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp
index 2d83d1091..c85c9c8d7 100644
--- a/src/Init/DiffusionInit.cpp
+++ b/src/Init/DiffusionInit.cpp
@@ -5,6 +5,7 @@
 #include <map>
 #include <set>
 
+#include "DataStructures/Field.hpp"
 #include "InitialList.hpp"
 
 #include <Rcpp/proxy/ProtectedProxy.h>
@@ -14,6 +15,7 @@
 #include <cstdint>
 #include <stdexcept>
 #include <string>
+#include <utility>
 #include <vector>
 
 namespace poet {
@@ -175,6 +177,42 @@ void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
   // R.parseEval("print(alpha_y)");
 }
 
+InitialList::DiffusionInit::BoundaryMap
+RcppListToBoundaryMap(const std::vector<std::string> &trans_names,
+                      const Rcpp::List &bound_list, uint32_t n_cols,
+                      uint32_t n_rows) {
+  InitialList::DiffusionInit::BoundaryMap map;
+
+  for (const auto &name : trans_names) {
+    const Rcpp::List &conc_list = bound_list[name];
+
+    tug::Boundary<TugType> bc(n_rows, n_cols);
+
+    for (const auto &[tug_index, r_init_name] : tug_side_mapping) {
+      const Rcpp::List &side_list = conc_list[tug_index];
+
+      const Rcpp::NumericVector &type = side_list["type"];
+      const Rcpp::NumericVector &value = side_list["value"];
+
+      if (type.size() != value.size()) {
+        throw std::runtime_error(
+            "Boundary type and value are not the same length");
+      }
+
+      for (std::size_t i = 0; i < type.size(); i++) {
+        if (type[i] == tug::BC_TYPE_CONSTANT) {
+          bc.setBoundaryElementConstant(static_cast<tug::BC_SIDE>(tug_index), i,
+                                        value[i]);
+        }
+      }
+    }
+
+    map[name] = bc.serialize();
+  }
+
+  return map;
+}
+
 InitialList::DiffusionInit InitialList::getDiffusionInit() const {
   DiffusionInit diff_init;
 
@@ -189,8 +227,8 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
   diff_init.constant_cells = this->constant_cells;
   diff_init.transport_names = this->transport_names;
 
-  diff_init.initial_grid = Field(this->initial_grid);
-  diff_init.boundaries = Field(this->boundaries);
+  diff_init.boundaries = RcppListToBoundaryMap(
+      this->transport_names, this->boundaries, this->n_cols, this->n_rows);
   diff_init.alpha_x = Field(this->alpha_x);
   diff_init.alpha_y = Field(this->alpha_y);
 
diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp
index 0931c34f3..63e3e2742 100644
--- a/src/Init/InitialList.cpp
+++ b/src/Init/InitialList.cpp
@@ -38,10 +38,10 @@ void InitialList::importList(const Rcpp::List &setup) {
       setup[static_cast<int>(ExportList::DIFFU_TRANSPORT)]);
   this->boundaries =
       Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_BOUNDARIES)]);
-  this->alpha_x = Rcpp::as<std::vector<double>>(
-      setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
-  this->alpha_y = Rcpp::as<std::vector<double>>(
-      setup[static_cast<int>(ExportList::DIFFU_ALPHA_Y)]);
+  this->alpha_x =
+      Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
+  this->alpha_y =
+      Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_Y)]);
 
   this->database =
       Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index 0f3b631b1..86c4ded77 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <tug/Boundary.hpp>
+
 #include <RInside.h>
 #include <Rcpp.h>
 #include <Rcpp/vector/instantiation.h>
@@ -8,7 +10,9 @@
 #include <cstdint>
 
 #include <IPhreeqcPOET.hpp>
+#include <map>
 #include <string>
+#include <utility>
 #include <vector>
 
 #include <DataStructures/Field.hpp>
@@ -26,6 +30,8 @@ public:
   void importList(const Rcpp::List &setup);
   Rcpp::List exportList();
 
+  Field getInitialGrid() const { return Field(this->initial_grid); }
+
 private:
   RInside &R;
 
@@ -93,6 +99,9 @@ private:
 
 public:
   struct DiffusionInit {
+    using BoundaryMap =
+        std::map<std::string, std::vector<tug::BoundaryElement<TugType>>>;
+
     uint8_t dim;
     std::uint32_t n_cols;
     std::uint32_t n_rows;
@@ -104,8 +113,8 @@ public:
 
     std::vector<std::string> transport_names;
 
-    Field initial_grid;
-    Field boundaries;
+    BoundaryMap boundaries;
+
     Field alpha_x;
     Field alpha_y;
   };
diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp
index abbecce95..17f761b57 100644
--- a/src/Transport/DiffusionModule.cpp
+++ b/src/Transport/DiffusionModule.cpp
@@ -4,6 +4,9 @@
 **
 ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
 **
+** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
+** of Potsdam)
+**
 ** POET is free software; you can redistribute it and/or modify it under the
 ** terms of the GNU General Public License as published by the Free Software
 ** Foundation; either version 2 of the License, or (at your option) any later
@@ -19,174 +22,101 @@
 */
 
 #include "DiffusionModule.hpp"
-#include "../Base/Macros.hpp"
 
-#include <tug/BoundaryCondition.hpp>
-#include <tug/Diffusion.hpp>
+#include "Base/Macros.hpp"
+#include "Init/InitialList.hpp"
 
-#include <Rcpp.h>
-#include <algorithm>
-#include <cstdint>
+#include <Eigen/src/Core/Map.h>
+#include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
+#include <chrono>
 
-#include <array>
-#include <cassert>
-#include <mpi.h>
+#include "tug/Boundary.hpp"
+#include "tug/Grid.hpp"
+#include "tug/Simulation.hpp"
+
+#include <iostream>
+#include <ostream>
 #include <string>
 #include <vector>
 
 using namespace poet;
 
-static constexpr double ZERO_MULTIPLICATOR = 10E-14;
+static inline std::vector<TugType>
+MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
+  std::vector<TugType> vec(mat.rows() * mat.cols());
 
-constexpr std::array<uint8_t, 4> borders = {
-    tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP,
-    tug::bc::BC_SIDE_BOTTOM};
-
-inline const char *convert_bc_to_config_file(uint8_t in) {
-  switch (in) {
-  case tug::bc::BC_SIDE_TOP:
-    return "N";
-  case tug::bc::BC_SIDE_RIGHT:
-    return "E";
-  case tug::bc::BC_SIDE_BOTTOM:
-    return "S";
-  case tug::bc::BC_SIDE_LEFT:
-    return "W";
-  }
-  return "";
-}
-
-DiffusionModule::DiffusionModule(const poet::DiffusionParams &diffu_args,
-                                 const poet::GridParams &grid_params)
-    : n_cells_per_prop(grid_params.total_n) {
-  this->diff_input.setGridCellN(grid_params.n_cells[0], grid_params.n_cells[1]);
-  this->diff_input.setDomainSize(grid_params.s_cells[0],
-                                 grid_params.s_cells[1]);
-
-  this->dim = grid_params.dim;
-
-  this->initialize(diffu_args, grid_params.total_n);
-}
-
-void DiffusionModule::initialize(const poet::DiffusionParams &args,
-                                 std::uint32_t n_grid_cells) {
-  // const poet::DiffusionParams args(this->R);
-
-  // name of props
-  //
-  this->prop_names = Rcpp::as<std::vector<std::string>>(args.initial_t.names());
-  this->prop_count = this->prop_names.size();
-
-  // initialize field and alphas
-  this->alpha.reserve(this->prop_count);
-
-  std::vector<std::vector<double>> initial_values;
-
-  for (uint32_t i = 0; i < this->prop_count; i++) {
-    // get alphas - we cannot assume alphas are provided in same order as
-    // initial input
-    this->alpha.push_back(args.alpha[this->prop_names[i]]);
-
-    const double val = args.initial_t[prop_names[i]];
-    std::vector<double> init_val(n_grid_cells, val);
-    initial_values.push_back(std::move(init_val));
-
-    if (this->dim == this->DIM_2D) {
-      tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0],
-                                    diff_input.grid.grid_cells[1]);
-      this->bc_vec.push_back(bc);
-    } else {
-      tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0]);
-      this->bc_vec.push_back(bc);
+  for (std::uint32_t i = 0; i < mat.cols(); i++) {
+    for (std::uint32_t j = 0; j < mat.rows(); j++) {
+      vec[j * mat.cols() + i] = mat(j, i);
     }
   }
 
-  this->t_field = Field(n_grid_cells, initial_values, prop_names);
-
-  // apply boundary conditions to each ghost node
-  uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4);
-  for (uint8_t i = 0; i < bc_size; i++) {
-    const auto &side = borders[i];
-    std::vector<uint32_t> vecinj_i = Rcpp::as<std::vector<uint32_t>>(
-        args.vecinj_index[convert_bc_to_config_file(side)]);
-    for (int i = 0; i < this->prop_count; i++) {
-      std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
-      tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
-      for (int j = 0; j < vecinj_i.size(); j++) {
-        if (vecinj_i[j] == 0) {
-          continue;
-        }
-        if (this->dim == this->DIM_2D) {
-          curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT,
-                              bc_vec[vecinj_i[j] - 1]};
-        }
-        if (this->dim == this->DIM_1D) {
-          curr_bc(side) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]};
-        }
-      }
-    }
-  }
-
-  // apply inner grid constant cells
-  // NOTE: opening a scope here for distinguish variable names
-  if (args.vecinj_inner.size() != 0) {
-    // apply inner constant cells for every concentration
-    for (int i = 0; i < this->prop_count; i++) {
-      std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
-      tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
-      for (int j = 0; j < args.vecinj_inner.size(); j++) {
-        std::vector<double> inner_tuple =
-            Rcpp::as<std::vector<double>>(args.vecinj_inner[j]);
-        tug::bc::boundary_condition bc = {tug::bc::BC_TYPE_CONSTANT,
-                                          bc_vec[inner_tuple[0] - 1]};
-
-        this->index_constant_cells.push_back(inner_tuple[1]);
-        uint32_t x = inner_tuple[1];
-        uint32_t y = (this->dim == this->DIM_1D ? 0 : inner_tuple[2]);
-
-        curr_bc.setInnerBC(bc, x, y);
-      }
-    }
-  }
+  return vec;
 }
 
-void DiffusionModule::simulate(double dt) {
-  double sim_a_transport, sim_b_transport;
+static inline Eigen::MatrixX<TugType>
+VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
+            std::uint32_t n_cols) {
+  Eigen::MatrixX<TugType> mat(n_rows, n_cols);
 
-  sim_b_transport = MPI_Wtime();
-
-  MSG_NOENDL("DiffusionModule::simulate(): Starting diffusion ...");
-  std::cout << std::flush;
-
-  std::vector<std::vector<double>> field_2d = t_field.As2DVector();
-
-  this->diff_input.setTimestep(dt);
-
-  for (int i = 0; i < field_2d.size(); i++) {
-    std::vector<double> in_alpha(this->n_cells_per_prop, this->alpha[i]);
-    this->diff_input.setBoundaryCondition(this->bc_vec[i]);
-
-    if (this->dim == this->DIM_1D) {
-      tug::diffusion::BTCS_1D(this->diff_input, field_2d[i].data(),
-                              in_alpha.data());
-    } else {
-      tug::diffusion::ADI_2D(this->diff_input, field_2d[i].data(),
-                             in_alpha.data());
+  for (std::uint32_t i = 0; i < n_cols; i++) {
+    for (std::uint32_t j = 0; j < n_rows; j++) {
+      mat(j, i) = vec[j * n_cols + i];
     }
   }
 
-  t_field = field_2d;
-
-  sim_a_transport = MPI_Wtime();
-
-  transport_t += sim_a_transport - sim_b_transport;
-  std::cout << " done in " << sim_a_transport - sim_b_transport << "sec"
-            << std::endl;
+  return mat;
 }
 
-void DiffusionModule::end() {
-  // R["simtime_transport"] = transport_t;
-  // R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
+// static constexpr double ZERO_MULTIPLICATOR = 10E-14;
+
+void DiffusionModule::simulate(double requested_dt) {
+  MSG("Starting diffusion ...");
+  const auto start_diffusion_t = std::chrono::high_resolution_clock::now();
+
+  const auto &n_rows = this->param_list.n_rows;
+  const auto &n_cols = this->param_list.n_cols;
+
+  tug::Grid<TugType> grid(param_list.n_rows, param_list.n_cols);
+  tug::Boundary boundary(grid);
+
+  grid.setDomain(param_list.s_rows, param_list.s_cols);
+
+  tug::Simulation<TugType> sim(grid, boundary);
+  sim.setIterations(1);
+
+  for (const auto &sol_name : this->param_list.transport_names) {
+    auto &species_conc = this->transport_field[sol_name];
+
+    Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
+    Eigen::MatrixX<TugType> alpha_x =
+        VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
+    Eigen::MatrixX<TugType> alpha_y =
+        VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols);
+
+    boundary.deserialize(this->param_list.boundaries[sol_name]);
+
+    grid.setAlpha(alpha_x, alpha_y);
+    grid.setConcentrations(conc);
+
+    sim.setTimestep(requested_dt);
+
+    sim.run();
+
+    species_conc = MatrixToVec(grid.getConcentrations());
+  }
+
+  const auto end_diffusion_t = std::chrono::high_resolution_clock::now();
+
+  const auto consumed_time_seconds =
+      std::chrono::duration_cast<std::chrono::duration<double>>(
+          end_diffusion_t - start_diffusion_t);
+
+  MSG("Diffusion done in " + std::to_string(consumed_time_seconds.count()) +
+      " [s]");
+
+  transport_t += consumed_time_seconds.count();
 }
 
 double DiffusionModule::getTransportTime() { return this->transport_t; }
diff --git a/src/Transport/DiffusionModule.hpp b/src/Transport/DiffusionModule.hpp
index a31bc0738..f19ee2064 100644
--- a/src/Transport/DiffusionModule.hpp
+++ b/src/Transport/DiffusionModule.hpp
@@ -21,20 +21,13 @@
 #ifndef DIFFUSION_MODULE_H
 #define DIFFUSION_MODULE_H
 
-#include "../Base/Grid.hpp"
-#include "../Base/SimParams.hpp"
-#include "../DataStructures/DataStructures.hpp"
+#include "DataStructures/Field.hpp"
+#include "Init/InitialList.hpp"
 
-#include <tug/BoundaryCondition.hpp>
-#include <tug/Diffusion.hpp>
-
-#include <array>
-#include <cmath>
-#include <cstdint>
-#include <string>
-#include <vector>
+#include <sys/types.h>
 
 namespace poet {
+
 /**
  * @brief Class describing transport simulation
  *
@@ -53,8 +46,8 @@ public:
    *
    * @param R RRuntime object
    */
-  DiffusionModule(const poet::DiffusionParams &diffu_args,
-                  const poet::GridParams &grid_params);
+  DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
+      : param_list(init_list), transport_field(field){};
 
   /**
    * @brief Run simulation for one iteration
@@ -64,14 +57,6 @@ public:
    */
   void simulate(double dt);
 
-  /**
-   * @brief End simulation
-   *
-   * All measured timings are distributed to the R runtime
-   *
-   */
-  void end();
-
   /**
    * @brief Get the transport time
    *
@@ -84,33 +69,17 @@ public:
    *
    * \return Reference to the diffusion field.
    */
-  Field &getField() { return this->t_field; }
+  Field &getField() { return this->transport_field; }
 
 private:
   /**
    * @brief Instance of RRuntime
    *
    */
-  // RRuntime &R;
 
-  enum { DIM_1D = 1, DIM_2D };
+  InitialList::DiffusionInit param_list;
 
-  void initialize(const poet::DiffusionParams &args,
-                  std::uint32_t n_grid_cells);
-
-  uint8_t dim;
-
-  uint32_t prop_count;
-
-  tug::diffusion::TugInput diff_input;
-  std::vector<double> alpha;
-  std::vector<uint32_t> index_constant_cells;
-  std::vector<std::string> prop_names;
-
-  std::vector<tug::bc::BoundaryCondition> bc_vec;
-  Field t_field;
-
-  uint32_t n_cells_per_prop;
+  Field transport_field;
 
   /**
    * @brief time spent for transport
diff --git a/src/poet.cpp b/src/poet.cpp
index 10c706a0d..fa7819460 100644
--- a/src/poet.cpp
+++ b/src/poet.cpp
@@ -4,6 +4,8 @@
 **
 ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
 **
+** Copyright (C) 2023-2024 Max Luebke (University of Potsdam)
+**
 ** POET is free software; you can redistribute it and/or modify it under the
 ** terms of the GNU General Public License as published by the Free Software
 ** Foundation; either version 2 of the License, or (at your option) any later
@@ -18,125 +20,308 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
-#include "Base/Grid.hpp"
 #include "Base/Macros.hpp"
 #include "Base/RInsidePOET.hpp"
-#include "Base/SimParams.hpp"
-#include "Chemistry/ChemistryModule.hpp"
+// #include "Chemistry/ChemistryModule.hpp"
+#include "DataStructures/Field.hpp"
+#include "Init/InitialList.hpp"
 #include "Transport/DiffusionModule.hpp"
 
+#include <RInside.h>
+#include <R_ext/Boolean.h>
+#include <Rcpp/DataFrame.h>
+#include <Rcpp/vector/instantiation.h>
 #include <poet.hpp>
 
 #include <Rcpp.h>
-#include <cstdint>
 #include <cstdlib>
 
-#include <cstring>
-#include <iostream>
 #include <string>
-#include <vector>
 
 #include <mpi.h>
 
+#include "Base/argh.hpp"
+
 using namespace std;
 using namespace poet;
 using namespace Rcpp;
 
-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());
+static int MY_RANK = 0;
 
-  for (const auto &name : col_names) {
-    double val = df[name.c_str()];
-    out_map.insert({name, val});
-  }
+// 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());
 
-  return out_map;
-}
+//   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) {
-  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)"));
+void writeFieldsToR(RInside &R, const Field &trans) {
+  // , const Field &chem) {
 
-  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)"));
+  Rcpp::DataFrame t_field(trans.asSEXP());
+
+  R["TMP"] = t_field;
+
+  R.parseEval("mysetup$state_T <- TMP");
+  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)"));
 }
 
-void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
-                         const std::string &database_path) {
-  chem.SetErrorHandlerMode(1);
-  chem.SetComponentH2O(false);
-  chem.SetRebalanceFraction(0.5);
-  chem.SetRebalanceByCell(true);
-  chem.UseSolutionDensityVolume(false);
-  chem.SetPartitionUZSolids(false);
+enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
 
-  // Set concentration units
-  // 1, mg/L; 2, mol/L; 3, kg/kgs
-  chem.SetUnitsSolution(2);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsPPassemblage(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsExchange(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsSurface(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsGasPhase(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsSSassemblage(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsKinetics(1);
+ParseRet parseInitValues(char **argv, RInsidePOET &R,
+                         RuntimeParameters &params) {
+  // initialize argh object
+  argh::parser cmdl(argv);
 
-  // Set representative volume
-  std::vector<double> rv;
-  rv.resize(wp_size, 1.0);
-  chem.SetRepresentativeVolume(rv);
+  // if user asked for help
+  if (cmdl[{"help", "h"}]) {
+    if (MY_RANK == 0) {
+      MSG("Todo");
+      MSG("See README.md for further information.");
+    }
 
-  // Set initial porosity
-  std::vector<double> por;
-  por.resize(wp_size, 1);
-  chem.SetPorosity(por);
+    return ParseRet::PARSER_HELP;
+  }
+  // if positional arguments are missing
+  else if (!cmdl(2)) {
+    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");
+    }
+    return ParseRet::PARSER_ERROR;
+  }
 
-  // Set initial saturation
-  std::vector<double> sat;
-  sat.resize(wp_size, 1.0);
-  chem.SetSaturation(sat);
+  // parse flags and check for unknown
+  for (const auto &option : cmdl.flags()) {
+    if (!(flaglist.find(option) != flaglist.end())) {
+      if (MY_RANK == 0) {
+        ERRMSG("Unrecognized option: " + option);
+        ERRMSG("Make sure to use available options. Exiting!");
+      }
+      return ParseRet::PARSER_ERROR;
+    }
+  }
 
-  // Load database
-  chem.LoadDatabase(database_path);
+  // parse parameters and check for unknown
+  for (const auto &option : cmdl.params()) {
+    if (!(paramlist.find(option.first) != paramlist.end())) {
+      if (MY_RANK == 0) {
+        ERRMSG("Unrecognized option: " + option.first);
+        ERRMSG("Make sure to use available options. Exiting!");
+      }
+      return ParseRet::PARSER_ERROR;
+    }
+  }
+
+  // simparams.print_progressbar = cmdl[{"P", "progress"}];
+
+  // // simparams.print_progressbar = cmdl[{"P", "progress"}];
+
+  // /* Parse DHT arguments */
+  // chem_params.use_dht = cmdl["dht"];
+  // chem_params.use_interp = cmdl["interp"];
+  // // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
+
+  // cmdl("dht-size", DHT_SIZE_PER_PROCESS_MB) >> chem_params.dht_size;
+  // // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
+  // // endl;
+
+  // cmdl("dht-snaps", 0) >> chem_params.dht_snaps;
+
+  // cmdl("dht-file") >> chem_params.dht_file;
+  // /*Parse work package size*/
+  // cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size;
+
+  // cmdl("interp-size", 100) >> chem_params.pht_size;
+  // cmdl("interp-min", 5) >> chem_params.interp_min_entries;
+  // cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
+
+  // /*Parse output options*/
+  // simparams.store_result = !cmdl["ignore-result"];
+
+  // /*Parse work package size*/
+  // cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size;
+
+  // chem_params.use_interp = cmdl["interp"];
+  // cmdl("interp-size", 100) >> chem_params.pht_size;
+  // cmdl("interp-min", 5) >> chem_params.interp_min_entries;
+  // cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
+
+  // /*Parse output options*/
+  // simparams.store_result = !cmdl["ignore-result"];
+
+  // if (simparams.world_rank == 0) {
+  //   MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
+  //   MSG("Work Package Size: " + std::to_string(simparams.wp_size));
+  //   MSG("DHT is " + BOOL_PRINT(chem_params.use_dht));
+
+  //   if (chem_params.use_dht) {
+  //     MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
+  //     // MDL: these should be outdated (?)
+  //     // MSG("DHT key default digits (ignored if 'signif_vector' is "
+  //     // 	"defined) = "
+  //     // 	 << simparams.dht_significant_digits);
+  //     // MSG("DHT logarithm before rounding: "
+  //     // 	 << (simparams.dht_log ? "ON" : "OFF"));
+  //     MSG("DHT size per process (Megabyte) = " +
+  //         std::to_string(chem_params.dht_size));
+  //     MSG("DHT save snapshots is " + BOOL_PRINT(chem_params.dht_snaps));
+  //     MSG("DHT load file is " + chem_params.dht_file);
+  //   }
+
+  //   if (chem_params.use_interp) {
+  //     MSG("PHT interpolation enabled: " +
+  //     BOOL_PRINT(chem_params.use_interp)); MSG("PHT interp-size = " +
+  //     std::to_string(chem_params.pht_size)); MSG("PHT interp-min  = " +
+  //         std::to_string(chem_params.interp_min_entries));
+  //     MSG("PHT interp-bucket-entries = " +
+  //         std::to_string(chem_params.pht_max_entries));
+  //   }
+  // }
+
+  std::string init_file;
+  std::string runtime_file;
+  std::string out_dir;
+
+  cmdl(1) >> init_file;
+  cmdl(2) >> runtime_file;
+  cmdl(3) >> out_dir;
+
+  // chem_params.dht_outdir = out_dir;
+
+  /* distribute information to R runtime */
+  // if local_rank == 0 then master else worker
+  R["local_rank"] = MY_RANK;
+  // assign a char* (string) to 'filesim'
+  R["filesim"] = wrap(runtime_file);
+  // assign a char* (string) to 'fileout'
+  R["fileout"] = wrap(out_dir);
+  // pass the boolean "store_result" to the R process
+  // R["store_result"] = simparams.store_result;
+  // // worker count
+  // R["n_procs"] = simparams.world_size - 1;
+  // // work package size
+  // R["work_package_size"] = simparams.wp_size;
+  // // dht enabled?
+  // R["dht_enabled"] = chem_params.use_dht;
+  // // log before rounding?
+  // R["dht_log"] = simparams.dht_log;
+
+  try {
+    Rcpp::Function source("source");
+    Rcpp::Function readRDS("readRDS");
+
+    Rcpp::List init_params_ = readRDS(init_file);
+    params.init_params = init_params_;
+
+    Rcpp::List runtime_params =
+        source(runtime_file, Rcpp::Named("local", true));
+    runtime_params = runtime_params["value"];
+    R[r_runtime_parameters] = runtime_params;
+
+    params.timesteps =
+        Rcpp::as<std::vector<double>>(runtime_params["timesteps"]);
+
+  } catch (const std::exception &e) {
+    ERRMSG("Error while parsing R scripts: " + std::string(e.what()));
+    return ParseRet::PARSER_ERROR;
+  }
+  // eval the init string, ignoring any returns
+  // R.parseEvalQ("source(filesim)");
+  // R.parseEvalQ("mysetup <- setup");
+
+  // this->chem_params.initFromR(R);
+
+  return ParseRet::PARSER_OK;
 }
 
-inline double RunMasterLoop(SimParams &params, RInside &R,
-                            const GridParams &g_params, uint32_t nxyz_master) {
+// void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
+//                          const std::string &database_path) {
+//   chem.SetErrorHandlerMode(1);
+//   chem.SetComponentH2O(false);
+//   chem.SetRebalanceFraction(0.5);
+//   chem.SetRebalanceByCell(true);
+//   chem.UseSolutionDensityVolume(false);
+//   chem.SetPartitionUZSolids(false);
+
+//   // Set concentration units
+//   // 1, mg/L; 2, mol/L; 3, kg/kgs
+//   chem.SetUnitsSolution(2);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsPPassemblage(1);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsExchange(1);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsSurface(1);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsGasPhase(1);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsSSassemblage(1);
+//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+//   chem.SetUnitsKinetics(1);
+
+//   // Set representative volume
+//   std::vector<double> rv;
+//   rv.resize(wp_size, 1.0);
+//   chem.SetRepresentativeVolume(rv);
+
+//   // Set initial porosity
+//   std::vector<double> por;
+//   por.resize(wp_size, 1);
+//   chem.SetPorosity(por);
+
+//   // Set initial saturation
+//   std::vector<double> sat;
+//   sat.resize(wp_size, 1.0);
+//   chem.SetSaturation(sat);
+
+//   // Load database
+//   chem.LoadDatabase(database_path);
+// }
+
+static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
+                            DiffusionModule &diffusion) {
 
-  DiffusionParams d_params{R};
-  DiffusionModule diffusion(d_params, g_params);
   /* Iteration Count is dynamic, retrieving value from R (is only needed by
    * master for the following loop) */
-  uint32_t maxiter = R.parseEval("mysetup$iterations");
+  uint32_t maxiter = params.timesteps.size();
 
   double sim_time = .0;
 
-  ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
-                       params.getChemParams(), MPI_COMM_WORLD);
+  // ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
+  //                      params.getChemParams(), MPI_COMM_WORLD);
 
-  set_chem_parameters(chem, nxyz_master, params.getChemParams().database_path);
-  chem.RunInitFile(params.getChemParams().input_script);
+  // set_chem_parameters(chem, nxyz_master,
+  // params.getChemParams().database_path);
+  // chem.RunInitFile(params.getChemParams().input_script);
 
-  poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t);
-  chem.initializeField(diffusion.getField());
+  // poet::ChemistryModule::SingleCMap init_df =
+  // DFToHashMap(d_params.initial_t);
+  // chem.initializeField(diffusion.getField());
 
-  if (params.getNumParams().print_progressbar) {
-    chem.setProgressBarPrintout(true);
-  }
+  // if (params.getNumParams().print_progressbar) {
+  //   chem.setProgressBarPrintout(true);
+  // }
 
   /* SIMULATION LOOP */
 
@@ -147,31 +332,30 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
     // cout << "CPP: Evaluating next time step" << endl;
     // R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
 
-    double dt = Rcpp::as<double>(
-        R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
+    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]");
 
     /* displaying iteration number, with C++ and R iterator */
     MSG("Going through iteration " + std::to_string(iter));
-    MSG("R's $iter: " +
-        std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
-        ". Iteration");
+    // MSG("R's $iter: " +
+    //     std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
+    //     ". Iteration");
 
     /* run transport */
     // TODO: transport to diffusion
     diffusion.simulate(dt);
 
-    chem.getField().update(diffusion.getField());
+    // chem.getField().update(diffusion.getField());
 
     MSG("Chemistry step");
 
-    chem.SetTimeStep(dt);
-    chem.RunCells();
+    // chem.SetTimeStep(dt);
+    // chem.RunCells();
 
-    writeFieldsToR(R, diffusion.getField(), chem.GetField());
-    diffusion.getField().update(chem.GetField());
+    writeFieldsToR(R, diffusion.getField());
+    // diffusion.getField().update(chem.GetField());
 
     R["req_dt"] = dt;
     R["simtime"] = (sim_time += dt);
@@ -183,7 +367,7 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
     // state_C after every iteration if the cmdline option
     // --ignore-results is not given (and thus the R variable
     // store_result is TRUE)
-    R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)");
+    R.parseEval("mysetup <- master_iteration_end(setup=mysetup)");
 
     MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
         std::to_string(maxiter));
@@ -194,60 +378,60 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
     dSimTime += end_t - start_t;
   } // END SIMULATION LOOP
 
-  R.parseEvalQ("profiling <- list()");
+  // R.parseEvalQ("profiling <- list()");
 
-  R["simtime_chemistry"] = chem.GetChemistryTime();
-  R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
+  // R["simtime_chemistry"] = chem.GetChemistryTime();
+  // R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
 
-  R["chemistry_loop"] = chem.GetMasterLoopTime();
-  R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop");
+  // R["chemistry_loop"] = chem.GetMasterLoopTime();
+  // R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop");
 
-  R["chemistry_sequential"] = chem.GetMasterSequentialTime();
-  R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential");
+  // R["chemistry_sequential"] = chem.GetMasterSequentialTime();
+  // R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential");
 
-  R["idle_master"] = chem.GetMasterIdleTime();
-  R.parseEvalQ("profiling$idle_master <- idle_master");
+  // R["idle_master"] = chem.GetMasterIdleTime();
+  // R.parseEvalQ("profiling$idle_master <- idle_master");
 
-  R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
-  R.parseEvalQ("profiling$idle_worker <- idle_worker");
+  // R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
+  // R.parseEvalQ("profiling$idle_worker <- idle_worker");
 
-  R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
-  R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
+  // R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
+  // R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
 
-  R["simtime_transport"] = diffusion.getTransportTime();
-  R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
+  // R["simtime_transport"] = diffusion.getTransportTime();
+  // R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
 
   // R["phreeqc_count"] = phreeqc_counts;
   // R.parseEvalQ("profiling$phreeqc_count <- phreeqc_count");
 
-  if (params.getChemParams().use_dht) {
-    R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
-    R.parseEvalQ("profiling$dht_hits <- dht_hits");
-    R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
-    R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
-    R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
-    R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
-    R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
-    R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
-  }
-  if (params.getChemParams().use_interp) {
-    R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
-    R.parseEvalQ("profiling$interp_write <- interp_w");
-    R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
-    R.parseEvalQ("profiling$interp_read <- interp_r");
-    R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
-    R.parseEvalQ("profiling$interp_gather <- interp_g");
-    R["interp_fc"] =
-        Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
-    R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
-    R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
-    R.parseEvalQ("profiling$interp_calls <- interp_calls");
-    R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
-    R.parseEvalQ("profiling$interp_cached <- interp_cached");
-  }
+  // if (params.getChemParams().use_dht) {
+  //   R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
+  //   R.parseEvalQ("profiling$dht_hits <- dht_hits");
+  //   R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
+  //   R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
+  //   R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
+  //   R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
+  //   R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
+  //   R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
+  // }
+  // if (params.getChemParams().use_interp) {
+  //   R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
+  //   R.parseEvalQ("profiling$interp_write <- interp_w");
+  //   R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
+  //   R.parseEvalQ("profiling$interp_read <- interp_r");
+  //   R["interp_g"] =
+  // Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
+  //   R.parseEvalQ("profiling$interp_gather <- interp_g");
+  //   R["interp_fc"] =
+  //       Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
+  //   R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
+  //   R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
+  //   R.parseEvalQ("profiling$interp_calls <- interp_calls");
+  //   R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
+  //   R.parseEvalQ("profiling$interp_cached <- interp_cached");
+  // }
 
-  chem.MasterLoopBreak();
-  diffusion.end();
+  // chem.MasterLoopBreak();
 
   return dSimTime;
 }
@@ -271,24 +455,24 @@ int main(int argc, char *argv[]) {
 
   if (world_rank > 0) {
     {
-      SimParams params(world_rank, world_size);
-      int pret = params.parseFromCmdl(argv, R);
+      // 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;
-      }
+      // 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);
+      // // 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();
+      // worker.WorkerLoop();
     }
 
     MPI_Barrier(MPI_COMM_WORLD);
@@ -304,33 +488,31 @@ int main(int argc, char *argv[]) {
   // TODO: kann raus
   R.parseEvalQ(kin_r_library);
 
-  SimParams params(world_rank, world_size);
-  int pret = params.parseFromCmdl(argv, R);
+  RuntimeParameters run_params;
 
-  if (pret == poet::PARSER_ERROR) {
+  switch (parseInitValues(argv, R, run_params)) {
+  case ParseRet::PARSER_ERROR:
+  case ParseRet::PARSER_HELP:
     MPI_Finalize();
-    return EXIT_FAILURE;
-  } else if (pret == poet::PARSER_HELP) {
-    MPI_Finalize();
-    return EXIT_SUCCESS;
+    return 0;
+  case ParseRet::PARSER_OK:
+    break;
   }
 
   MSG("RInside initialized on process " + std::to_string(world_rank));
 
-  R.parseEvalQ("mysetup <- setup");
-  // if (world_rank == 0) { // get timestep vector from
-  // grid_init function ... //
-  std::string master_init_code = "mysetup <- master_init(setup=setup)";
+  // R.parseEvalQ("mysetup <- setup");
+  // // if (world_rank == 0) { // get timestep vector from
+  // // grid_init function ... //
+  std::string master_init_code = "mysetup <- master_init(setup=mysetup)";
   R.parseEval(master_init_code);
 
-  GridParams g_params(R);
-
-  params.initVectorParams(R);
+  // run_params.initVectorParams(R);
 
   // MDL: store all parameters
   if (world_rank == 0) {
     MSG("Calling R Function to store calling parameters");
-    R.parseEvalQ("StoreSetup(setup=mysetup)");
+    // R.parseEvalQ("StoreSetup(setup=mysetup)");
   }
 
   if (world_rank == 0) {
@@ -339,23 +521,28 @@ int main(int argc, char *argv[]) {
 
   // MPI_Barrier(MPI_COMM_WORLD);
 
-  uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1);
+  InitialList init_list(R);
+  init_list.importList(run_params.init_params);
 
-  dSimTime = RunMasterLoop(params, R, g_params, nxyz_master);
+  DiffusionModule diffusion(init_list.getDiffusionInit(),
+                            init_list.getInitialGrid());
+
+  dSimTime = RunMasterLoop(R, run_params, diffusion);
 
   MSG("finished simulation loop");
 
   MSG("start timing profiling");
 
-  R["simtime"] = dSimTime;
-  R.parseEvalQ("profiling$simtime <- simtime");
+  // R["simtime"] = dSimTime;
+  // R.parseEvalQ("profiling$simtime <- simtime");
 
-  string r_vis_code;
-  r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
-  R.parseEval(r_vis_code);
+  // string r_vis_code;
+  // r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
+  // R.parseEval(r_vis_code);
 
-  MSG("Done! Results are stored as R objects into <" + params.getOutDir() +
-      "/timings.rds>");
+  // MSG("Done! Results are stored as R objects into <" + run_params.getOutDir()
+  // +
+  //     "/timings.rds>");
 
   MPI_Barrier(MPI_COMM_WORLD);
 
diff --git a/src/poet.hpp.in b/src/poet.hpp.in
index 115ac7dd7..99b3e7a58 100644
--- a/src/poet.hpp.in
+++ b/src/poet.hpp.in
@@ -1,12 +1,79 @@
-#ifndef POET_H
-#define POET_H
+/*
+** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
+** Potsdam)
+**
+** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
+**
+** Copyright (C) 2023-2024 Max Luebke (University of Potsdam)
+**
+** POET is free software; you can redistribute it and/or modify it under the
+** terms of the GNU General Public License as published by the Free Software
+** Foundation; either version 2 of the License, or (at your option) any later
+** version.
+**
+** POET is distributed in the hope that it will be useful, but WITHOUT ANY
+** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License along with
+** this program; if not, write to the Free Software Foundation, Inc., 51
+** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
 
+#pragma once
+
+#include <set>
 #include <string>
+#include <vector>
+
+#include <Rcpp.h>
+
 static const char *poet_version = "@POET_VERSION@";
 
 // using the Raw string literal to avoid escaping the quotes
-static inline std::string kin_r_library = R"(@R_KIN_LIB@)";
+static const inline std::string kin_r_library = R"(@R_KIN_LIB@)";
 
-static inline std::string init_r_library = R"(@R_INIT_LIB@)";
+static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
 
-#endif // POET_H
+static const inline std::string r_runtime_parameters = "mysetup";
+
+const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
+                                     "interp"};
+const std::set<std::string> paramlist{
+    "work-package-size", "dht-strategy", "dht-size",   "dht-snaps",
+    "dht-file",          "interp-size",  "interp-min", "interp-bucket-entries"};
+
+struct RuntimeParameters {
+  std::vector<double> timesteps;
+
+  bool print_progressbar;
+
+  Rcpp::List init_params;
+
+  struct ChemistryParams {
+    // std::string database_path;
+    // std::string input_script;
+
+    //   bool use_dht;
+    //   std::uint64_t dht_size;
+    //   int dht_snaps;
+    //   std::string dht_file;
+    //   std::string dht_outdir;
+    //   NamedVector<std::uint32_t> dht_signifs;
+
+    //   bool use_interp;
+    //   std::uint64_t pht_size;
+    //   std::uint32_t pht_max_entries;
+    //   std::uint32_t interp_min_entries;
+    //   NamedVector<std::uint32_t> pht_signifs;
+
+    //   struct Chem_Hook_Functions {
+    //     RHookFunction<bool> dht_fill;
+    //     RHookFunction<std::vector<double>> dht_fuzz;
+    //     RHookFunction<std::vector<std::size_t>> interp_pre;
+    //     RHookFunction<bool> interp_post;
+    //   } hooks;
+
+    //   void initFromR(RInsidePOET &R);
+  };
+};