diff --git a/CMake/POET_Scripts.cmake b/CMake/POET_Scripts.cmake
index 89fc1e84f..13b66c545 100644
--- a/CMake/POET_Scripts.cmake
+++ b/CMake/POET_Scripts.cmake
@@ -13,7 +13,7 @@ macro(get_POET_version)
       WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
       OUTPUT_VARIABLE POET_GIT_VERSION
       OUTPUT_STRIP_TRAILING_WHITESPACE)
-    if(NOT POET_GIT_BRANCH STREQUAL "master")
+    if(NOT POET_GIT_BRANCH STREQUAL "main")
       set(POET_VERSION "${POET_GIT_BRANCH}/${POET_GIT_VERSION}")
     else()
       set(POET_VERSION "${POET_GIT_VERSION}")
@@ -21,7 +21,7 @@ macro(get_POET_version)
   elseif(EXISTS ${PROJECT_SOURCE_DIR}/.svn)
     file(STRINGS .gitversion POET_VERSION)
   else()
-    set(POET_VERSION "0.1")
+    set(POET_VERSION "not_versioned")
   endif()
 
   message(STATUS "Configuring POET version ${POET_VERSION}")
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ce0382d45..fa7f009a1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -41,4 +41,4 @@ endif()
 
 option(BUILD_DOC "Build documentation with doxygen" OFF)
 
-add_subdirectory(docs)
\ No newline at end of file
+add_subdirectory(docs)
diff --git a/README.md b/README.md
index b17d73c12..598117675 100644
--- a/README.md
+++ b/README.md
@@ -1,53 +1,63 @@
-<!--
-    Time-stamp: "Last modified 2023-08-02 13:55:11 mluebke"
--->
-
 # POET
 
-[POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive transport
-simulator implementing a parallel architecture and a fast, original MPI-based
-Distributed Hash Table.
+[POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive
+transport simulator implementing a parallel architecture and a fast,
+original MPI-based Distributed Hash Table.
 
 ![POET's Coupling Scheme](./docs/Scheme_POET_en.svg)
 
 ## Parsed code documentiation
 
-A parsed version of POET's documentiation can be found at [Gitlab
+A parsed version of POET's documentation can be found at [Gitlab
 pages](https://naaice.git-pages.gfz-potsdam.de/poet).
 
 ## External Libraries
 
-The following external header library is shipped with POET:
+The following external libraries are shipped with POET:
 
-- **argh** - https://github.com/adishavit/argh (BSD license)
-- **IPhreeqc** with patches from GFZ -
-  https://github.com/usgs-coupled/iphreeqc -
-  https://git.gfz-potsdam.de/naaice/iphreeqc
-- **tug** - https://git.gfz-potsdam.de/naaice/tug
+- **CLI11** - <https://github.com/CLIUtils/CLI11>
+- **IPhreeqc** with patches from GFZ/UP -
+  <https://github.com/usgs-coupled/iphreeqc> -
+  <https://git.gfz-potsdam.de/naaice/iphreeqc>
+- **tug** - <https://git.gfz-potsdam.de/naaice/tug>
 
 ## Installation
 
 ### Requirements
 
-To compile POET you need several software to be installed:
+To compile POET you need following software to be installed:
 
 - C/C++ compiler (tested with GCC)
 - MPI-Implementation (tested with OpenMPI and MVAPICH)
-- R language and environment
 - CMake 3.9+
 - Eigen3 3.4+ (required by `tug`)
-- *optional*: `doxygen` with `dot` bindings for documentiation
+- *optional*: `doxygen` with `dot` bindings for documentation
+- R language and environment including headers or `-dev` packages
+  (distro dependent)
 
-The following R libraries must then be installed, which will get the
-needed dependencies automatically:
+The following R packages (and their dependencies) must also be
+installed:
 
 - [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html)
 - [RInside](https://cran.r-project.org/web/packages/RInside/index.html)
+- [qs](https://cran.r-project.org/web/packages/qs/index.html)
+
+This can be simply achieved by issuing the following commands:
+
+```sh
+# start R environment
+$ R
+
+# install R dependencies (case sensitive!)
+> install.packages(c("Rcpp", "RInside","qs"))
+> q(save="no")
+```
+
 
 ### Compiling source code
 
-The generation of makefiles is done with CMake. You should be able to generate
-Makefiles by running:
+POET is built with CMake. You can generate Makefiles by running the
+usual:
 
 ```sh
 mkdir build && cd build
@@ -58,23 +68,23 @@ This will create the directory `build` and processes the CMake files
 and generate Makefiles from it. You're now able to run `make` to start
 build process.
 
-If everything went well you'll find the executable at
-`build/app/poet`, but it is recommended to install the POET project
+If everything went well you'll find the executables at
+`build/src/poet`, but it is recommended to install the POET project
 structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
 
 During the generation of Makefiles, various options can be specified
 via `cmake -D <option>=<value> [...]`. Currently, there are the
 following available options:
 
-- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed statistics about
-  DHT usage. Defaults to _OFF_.
-- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests (more to
-  come). Defaults to _OFF_.
-- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of accesses to one
-  PHT bucket. Use with caution, as things will get slowed down significantly.
-  Defaults to _OFF_.
-- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of predefined
-  models/benchmarks. Defaults to *ON*.
+- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed
+  statistics about DHT usage. Defaults to _OFF_.
+- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests
+  (more to come). Defaults to _OFF_.
+- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of
+  accesses to one PHT bucket. Use with caution, as things will get
+  slowed down significantly. Defaults to _OFF_.
+- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of
+  predefined models/benchmarks. Defaults to *ON*.
   
 ### Example: Build from scratch
 
@@ -87,7 +97,7 @@ follows:
 $ R
 
 # install R dependencies
-> install.packages(c("Rcpp", "RInside"))
+> install.packages(c("Rcpp", "RInside","qs"))
 > q(save="no")
 
 # cd into POET project root
@@ -131,38 +141,42 @@ poet
 
 With the installation of POET, two executables are provided: 
   - `poet` - the main executable to run simulations
-  - `poet_init` - a preprocessor to generate input files for POET from R scripts
+  - `poet_init` - a preprocessor to generate input files for POET from
+    R scripts
 
-Preprocessed benchmarks can be found in the `share/poet` directory with an
-according *runtime* setup. More on those files and how to create them later. 
+Preprocessed benchmarks can be found in the `share/poet` directory
+with an according *runtime* setup. More on those files and how to
+create them later.
 
 ## Running
 
-Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE> <OUTPUT_DIRECTORY>`
-where:
+Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE>
+<OUTPUT_DIRECTORY>` where:
 
 - **OPTIONS** - POET options (explained below)
-- **RUNFILE** - Runtime parameters described as R script 
+- **RUNFILE** - Runtime parameters described as R script
 - **SIMFILE** - Simulation input prepared by `poet_init`
-- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
+- **OUTPUT_DIRECTORY** - path, where all output of POET should be
+  stored
 
-### POET options
+### POET command line arguments
 
 The following parameters can be set:
 
-| Option                      | Value        | Description                                                                                                              |
-|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
-| **--work-package-size=**    | _1..n_       | size of work packages (defaults to _5_)                                                                                  |
-| **-P, --progress**          |              | show progress bar                                                                                                        |
-| **--ai-surrogate**          |              | activates the AI surrogate chemistry model (defaults to _OFF_) |
-| **--dht**                   |              | enabling DHT usage (defaults to _OFF_)                                                                                   |
-| **--dht-strategy=**         | _0-1_        | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_)                                                           |
-| **--dht-size=**             | _1-n_        | size of DHT per process involved in megabyte (defaults to _1000 MByte_)                                                  |
-| **--dht-snaps=**            | _0-2_        | disable or enable storage of DHT snapshots                                                                               |
-| **--dht-file=**             | `<SNAPSHOT>` | initializes DHT with the given snapshot file                                                                             |
-| **--interp-size**           | _1-n_        | size of PHT (interpolation) per process in megabyte                                                                      |
-| **--interp-bucket-entries** | _1-n_        | number of entries to store at maximum in one PHT bucket                                                                  |
-| **--interp-min**            | _1-n_        | number of entries in PHT bucket needed to start interpolation                                                            |
+| Option                      | Value        | Description                                                                      |
+|-----------------------------|--------------|----------------------------------------------------------------------------------|
+| **--work-package-size=**    | _1..n_       | size of work packages (defaults to _5_)                                          |
+| **-P, --progress**          |              | show progress bar                                                                |
+| **--ai-surrogate**          |              | activates the AI surrogate chemistry model (defaults to _OFF_)                   |
+| **--dht**                   |              | enabling DHT usage (defaults to _OFF_)                                           |
+| **--qs**                    |              | store results using qs::qsave() (.qs extension) instead of default RDS (.rds)    |
+| **--dht-strategy=**         | _0-1_        | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_)                   |
+| **--dht-size=**             | _1-n_        | size of DHT per process involved in megabyte (defaults to _1000 MByte_)          |
+| **--dht-snaps=**            | _0-2_        | disable or enable storage of DHT snapshots                                       |
+| **--dht-file=**             | `<SNAPSHOT>` | initializes DHT with the given snapshot file                                     |
+| **--interp-size**           | _1-n_        | size of PHT (interpolation) per process in megabyte                              |
+| **--interp-bucket-entries** | _1-n_        | number of entries to store at maximum in one PHT bucket                          |
+| **--interp-min**            | _1-n_        | number of entries in PHT bucket needed to start interpolation                    |
 
 #### Additions to `dht-snaps`
 
@@ -176,12 +190,13 @@ Following values can be set:
 
 ### Example: Running from scratch
 
-We will continue the above example and start a simulation with *barite_het*,
-which simulation files can be found in
-`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a heterogeneous
-diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with
-a time step size of 100 seconds. To start the simulation with 4 processes `cd`
-into your previously installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
+We will continue the above example and start a simulation with
+*barite_het*, which simulation files can be found in
+`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a
+heterogeneous diffusion is used. It's a small 2D grid, 2x5 grid,
+simulating 50 time steps with a time step size of 100 seconds. To
+start the simulation with 4 processes `cd` into your previously
+installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
 
 ```sh
 cp ../share/poet/barite/barite_het* .
@@ -191,11 +206,11 @@ mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output
 After a finished simulation all data generated by POET will be found
 in the directory `output`.
 
-You might want to use the DHT to cache previously simulated data and reuse them
-in further time-steps. Just append `--dht` to the options of POET to activate
-the usage of the DHT. Also, after each iteration a DHT snapshot shall be
-produced. This is done by appending the `--dht-snaps=<value>` option. The
-resulting call would look like this:
+You might want to use the DHT to cache previously simulated data and
+reuse them in further time-steps. Just append `--dht` to the options
+of POET to activate the usage of the DHT. Also, after each iteration a
+DHT snapshot shall be produced. This is done by appending the
+`--dht-snaps=<value>` option. The resulting call would look like this:
 
 ```sh
 mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.rds output
@@ -253,12 +268,13 @@ produce any valid predictions.
 
 ## Defining a model
 
-In order to provide a model to POET, you need to setup a R script which can then
-be used by `poet_init` to generate the simulation input. Which parameters are
-required can be found in the
-[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization). We try to
-keep the document up-to-date. However, if you encounter missing information or
-need help, please get in touch with us via the issue tracker or E-Mail.
+In order to provide a model to POET, you need to setup a R script
+which can then be used by `poet_init` to generate the simulation
+input. Which parameters are required can be found in the
+[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization).
+We try to keep the document up-to-date. However, if you encounter
+missing information or need help, please get in touch with us via the
+issue tracker or E-Mail.
 
 `poet_init` can be used as follows:
 
@@ -268,46 +284,59 @@ need help, please get in touch with us via the issue tracker or E-Mail.
 
 where: 
 
-- **output** - name of the output file (defaults to the input file name
-  with the extension `.rds`)
-- **setwd** - set the working directory to the directory of the input file (e.g.
-  to allow relative paths in the input script). However, the output file
-  will be stored in the directory from which `poet_init` was called.
+- **output** - name of the output file (defaults to the input file
+  name with the extension `.rds`)
+- **setwd** - set the working directory to the directory of the input
+  file (e.g. to allow relative paths in the input script). However,
+  the output file will be stored in the directory from which
+  `poet_init` was called.
+
+## About the usage of MPI_Wtime()
+
+Implemented time measurement functions uses `MPI_Wtime()`. Some
+important information from the OpenMPI Man Page:
+
+For example, on platforms that support it, the clock_gettime()
+function will be used to obtain a monotonic clock value with whatever
+precision is supported on that platform (e.g., nanoseconds).
 
 ## Additional functions for the AI surrogate
 
-The AI surrogate can be activated for any benchmark and is by default initiated
-as a sequential keras model with three hidden layer of depth 48, 96, 24 with
-relu activation and adam optimizer. All functions in `ai_surrogate_model.R` can
-be overridden by adding custom definitions via an R file in the input script.
-This is done by adding the path to this file in the input script. Simply add the
-path as an element called `ai_surrogate_input_script` to the `chemistry_setup`
-list. Please use the global variable `ai_surrogate_base_path` as a base path
+The AI surrogate can be activated for any benchmark and is by default
+initiated as a sequential keras model with three hidden layer of depth
+48, 96, 24 with relu activation and adam optimizer. All functions in
+`ai_surrogate_model.R` can be overridden by adding custom definitions
+via an R file in the input script. This is done by adding the path to
+this file in the input script. Simply add the path as an element
+called `ai_surrogate_input_script` to the `chemistry_setup` list.
+Please use the global variable `ai_surrogate_base_path` as a base path
 when relative filepaths are used in custom funtions.
 
-**There is currently no default implementation to determine the validity of
-predicted values.** This means, that every input script must include an R source
-file with a custom function `validate_predictions(predictors, prediction)`.
-Examples for custom functions can be found for the barite_200 benchmark
+**There is currently no default implementation to determine the
+validity of predicted values.** This means, that every input script
+must include an R source file with a custom function
+`validate_predictions(predictors, prediction)`. Examples for custom
+functions can be found for the barite_200 benchmark
 
-The functions can be defined as follows:  
+The functions can be defined as follows:
 
-`validate_predictions(predictors, prediction)`: Returns a boolean index vector
-that signals for each row in the predictions if the values are considered valid.
-Can eg. be implemented as a mass balance threshold between the predictors and
-the prediction.
+`validate_predictions(predictors, prediction)`: Returns a boolean
+index vector that signals for each row in the predictions if the
+values are considered valid. Can eg. be implemented as a mass balance
+threshold between the predictors and the prediction.
 
-`initiate_model()`: Returns a keras model. Can be used to load pretrained
-models.
+`initiate_model()`: Returns a keras model. Can be used to load
+pretrained models.
 
 `preprocess(df, backtransform = FALSE, outputs = FALSE)`: Returns the
-scaled/transformed/backtransformed dataframe. The `backtransform` flag signals
-if the current processing step is applied to data that's assumed to be scaled
-and expects backtransformed values. The `outputs` flag signals if the current
-processing step is applied to the output or tatget of the model. This can be
-used to eg. skip these processing steps and only scale the model input.
+scaled/transformed/backtransformed dataframe. The `backtransform` flag
+signals if the current processing step is applied to data that's
+assumed to be scaled and expects backtransformed values. The `outputs`
+flag signals if the current processing step is applied to the output
+or tatget of the model. This can be used to eg. skip these processing
+steps and only scale the model input.
 
-`training_step (model, predictor, target, validity)`: Trains the model after
-each iteration. `validity` is the bool index vector given by
-`validate_predictions` and can eg. be used to only train on values that have not
-been valid predictions.
\ No newline at end of file
+`training_step (model, predictor, target, validity)`: Trains the model
+after each iteration. `validity` is the bool index vector given by
+`validate_predictions` and can eg. be used to only train on values
+that have not been valid predictions.
diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R
index d0971e9da..468b485bc 100644
--- a/R_lib/init_r_lib.R
+++ b/R_lib/init_r_lib.R
@@ -1,15 +1,33 @@
-pqc_to_grid <- function(pqc_in, grid) {
+### Copyright (C) 2018-2024 Marco De Lucia, Max Luebke (GFZ Potsdam, 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.
+
+##' @param pqc_mat matrix, containing IDs and PHREEQC outputs 
+##' @param grid matrix, zonation referring to pqc_mat$ID 
+##' @return a data.frame
+pqc_to_grid <- function(pqc_mat, grid) {
     # Convert the input DataFrame to a matrix
-    pqc_in <- as.matrix(pqc_in)
+    pqc_mat <- as.matrix(pqc_mat)
 
     # Flatten the matrix into a vector
-    id_vector <- as.numeric(t(grid))
+    id_vector <- as.integer(t(grid))
 
     # Find the matching rows in the matrix
-    row_indices <- match(id_vector, pqc_in[, "ID"])
+    row_indices <- match(id_vector, pqc_mat[, "ID"])
 
-    # Extract the matching rows from pqc_in to size of grid matrix
-    result_mat <- pqc_in[row_indices, ]
+    # Extract the matching rows from pqc_mat to size of grid matrix
+    result_mat <- pqc_mat[row_indices, ]
 
     # Convert the result matrix to a data frame
     res_df <- as.data.frame(result_mat)
@@ -23,6 +41,12 @@ pqc_to_grid <- function(pqc_in, grid) {
     return(res_df)
 }
 
+
+##' @param pqc_mat matrix, 
+##' @param transport_spec column name of species in pqc_mat
+##' @param id
+##' @title 
+##' @return 
 resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
     df <- as.data.frame(pqc_mat, check.names = FALSE)
     value <- df[df$ID == id, transport_spec]
@@ -34,6 +58,10 @@ resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
     return(value)
 }
 
+##' @title 
+##' @param init_grid 
+##' @param new_names 
+##' @return 
 add_missing_transport_species <- function(init_grid, new_names) {
     # add 'ID' to new_names front, as it is not a transport species but required
     new_names <- c("ID", new_names)
diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R
index cb8eaecd3..cb5552a9b 100644
--- a/R_lib/kin_r_library.R
+++ b/R_lib/kin_r_library.R
@@ -1,5 +1,3 @@
-## Time-stamp: "Last modified 2023-08-15 11:58:23 delucia"
-
 ### 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
@@ -35,14 +33,18 @@ master_init <- function(setup, out_dir, init_field) {
     setup$iterations <- setup$maxiter
     setup$simulation_time <- 0
 
+    dgts <- as.integer(ceiling(log10(setup$maxiter)))
+    ## string format to use in sprintf
+    fmt <- paste0("%0", dgts, "d")
+
     if (is.null(setup[["store_result"]])) {
         setup$store_result <- TRUE
     }
 
     if (setup$store_result) {
-        init_field_out <- paste0(out_dir, "/iter_0.rds")
+        init_field_out <- paste0(out_dir, "/iter_", sprintf(fmt = fmt, 0), ".", setup$out_ext)
         init_field <- data.frame(init_field, check.names = FALSE)
-        saveRDS(init_field, file = init_field_out)
+        SaveRObj(x = init_field, path = init_field_out)
         msgm("Stored initial field in ", init_field_out)
         if (is.null(setup[["out_save"]])) {
             setup$out_save <- seq(1, setup$iterations)
@@ -61,7 +63,7 @@ master_iteration_end <- function(setup, state_T, state_C) {
     iter <- setup$iter
     # print(iter)
     ## max digits for iterations
-    dgts <- as.integer(ceiling(log10(setup$maxiter)))
+    dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
     ## string format to use in sprintf
     fmt <- paste0("%0", dgts, "d")
 
@@ -69,21 +71,23 @@ master_iteration_end <- function(setup, state_T, state_C) {
     ## comprised in setup$out_save
     if (setup$store_result) {
         if (iter %in% setup$out_save) {
-            nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
+            nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".", setup$out_ext)
             state_T <- data.frame(state_T, check.names = FALSE)
             state_C <- data.frame(state_C, check.names = FALSE)
-            
+
             ai_surrogate_info <- list(
-                prediction_time = if(exists("ai_prediction_time")) as.integer(ai_prediction_time) else NULL,
-                training_time = if(exists("ai_training_time")) as.integer(ai_training_time) else NULL,
-                valid_predictions = if(exists("validity_vector")) validity_vector else NULL)
-            saveRDS(list(
+                prediction_time = if (exists("ai_prediction_time")) as.integer(ai_prediction_time) else NULL,
+                training_time = if (exists("ai_training_time")) as.integer(ai_training_time) else NULL,
+                valid_predictions = if (exists("validity_vector")) validity_vector else NULL
+            )
+
+            SaveRObj(x = list(
                 T = state_T,
                 C = state_C,
                 simtime = as.integer(setup$simulation_time),
                 totaltime = as.integer(totaltime),
                 ai_surrogate_info = ai_surrogate_info
-            ), file = nameout)
+            ), path = nameout)
             msgm("results stored in <", nameout, ">")
         }
     }
@@ -172,3 +176,31 @@ GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
     ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
     return(as.integer(table(ids)))
 }
+
+
+## Handler to read R objs from binary files using either builtin
+## readRDS() or qs::qread() based on file extension
+ReadRObj <- function(path) {
+    ## code borrowed from tools::file_ext()
+    pos <- regexpr("\\.([[:alnum:]]+)$", path)
+    extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
+
+    switch(extension,
+        rds = readRDS(path),
+        qs  = qs::qread(path)
+    )
+}
+
+## Handler to store R objs to binary files using either builtin
+## saveRDS() or qs::qsave() based on file extension
+SaveRObj <- function(x, path) {
+    msgm("Storing to", path)
+    ## code borrowed from tools::file_ext()
+    pos <- regexpr("\\.([[:alnum:]]+)$", path)
+    extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
+
+    switch(extension,
+        rds = saveRDS(object = x, file = path),
+        qs  = qs::qsave(x = x, file = path)
+    )
+}
diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt
index d407713ba..01dc43caf 100644
--- a/bench/CMakeLists.txt
+++ b/bench/CMakeLists.txt
@@ -6,18 +6,19 @@ function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH)
 
     foreach(BENCH_FILE ${${POET_BENCH_LIST}})
         get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE)
-        set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME}.rds)
+        set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME})
+        set(OUT_FILE_EXT ${OUT_FILE}.qs)
 
         add_custom_command(
-            OUTPUT ${OUT_FILE}
-            COMMAND $<TARGET_FILE:poet_init> -o ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
+            OUTPUT ${OUT_FILE_EXT}
+            COMMAND $<TARGET_FILE:poet_init> -n ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
             COMMENT "Running poet_init on ${BENCH_FILE}"
             DEPENDS poet_init ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
             VERBATIM
             COMMAND_EXPAND_LISTS
         )
 
-        list(APPEND OUT_FILES_LIST ${OUT_FILE})
+        list(APPEND OUT_FILES_LIST ${OUT_FILE_EXT})
 
     endforeach(BENCH_FILE ${${POET_BENCH_LIST}})
         
diff --git a/ext/iphreeqc b/ext/iphreeqc
index 48e65d87a..e6e5e0d51 160000
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@@ -1 +1 @@
-Subproject commit 48e65d87ad70f84aec01c27d9560cd3094a8129c
+Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26
diff --git a/src/Base/Macros.hpp b/src/Base/Macros.hpp
index 47beeeae3..cb62d9749 100644
--- a/src/Base/Macros.hpp
+++ b/src/Base/Macros.hpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-09 14:16:04 mluebke"
 #ifndef MACROS_H
 #define MACROS_H
 
diff --git a/src/Base/RInsidePOET.hpp b/src/Base/RInsidePOET.hpp
index 2897fc5a8..466c49375 100644
--- a/src/Base/RInsidePOET.hpp
+++ b/src/Base/RInsidePOET.hpp
@@ -1,17 +1,13 @@
-#ifndef RPOET_H_
-#define RPOET_H_
+#pragma once
 
 #include <RInside.h>
 #include <Rcpp.h>
 #include <Rinternals.h>
-#include <cstddef>
 #include <exception>
-#include <optional>
-#include <stdexcept>
+#include <memory>
 #include <string>
-#include <utility>
-#include <vector>
 
+namespace poet {
 class RInsidePOET : public RInside {
 public:
   static RInsidePOET &getInstance() {
@@ -33,44 +29,64 @@ private:
   RInsidePOET() : RInside(){};
 };
 
-template <typename T> class RHookFunction {
+/**
+ * @brief Deferred evaluation function
+ *
+ * The class is intended to call R functions within an existing RInside
+ * instance. The problem with "original" Rcpp::Function is that they require:
+ * 1. RInside instance already present, restricting the declaration of
+ *     Rcpp::Functions in global scope
+ * 2. Require the function to be present. Otherwise, they will throw an
+ *     exception.
+ * This class solves both problems by deferring the evaluation of the function
+ *  until the constructor is called and evaluating whether the function is
+ *  present or not, wihout throwing an exception.
+ *
+ * @tparam T Return type of the function
+ */
+class DEFunc {
 public:
-  RHookFunction() {}
-  RHookFunction(RInside &R, const std::string &f_name) {
+  DEFunc() {}
+  DEFunc(const std::string &f_name) {
     try {
-      this->func = Rcpp::Function(Rcpp::as<SEXP>(R.parseEval(f_name.c_str())));
+      this->func = std::make_shared<Rcpp::Function>(f_name);
     } catch (const std::exception &e) {
     }
   }
 
-  RHookFunction(SEXP f) {
+  DEFunc(SEXP f) {
     try {
-      this->func = Rcpp::Function(f);
+      this->func = std::make_shared<Rcpp::Function>(f);
     } catch (const std::exception &e) {
     }
   }
 
-  template <typename... Args> T operator()(Args... args) const {
-    if (func.has_value()) {
-      return (Rcpp::as<T>(this->func.value()(args...)));
+  template <typename... Args> SEXP operator()(Args... args) const {
+    if (func) {
+      return (*this->func)(args...);
     } else {
       throw std::exception();
     }
   }
 
-  RHookFunction &operator=(const RHookFunction &rhs) {
+  DEFunc &operator=(const DEFunc &rhs) {
     this->func = rhs.func;
     return *this;
   }
 
-  RHookFunction(const RHookFunction &rhs) { this->func = rhs.func; }
+  DEFunc(const DEFunc &rhs) { this->func = rhs.func; }
 
-  bool isValid() const { return this->func.has_value(); }
+  bool isValid() const { return static_cast<bool>(func); }
 
-  SEXP asSEXP() const { return Rcpp::as<SEXP>(this->func.value()); }
+  SEXP asSEXP() const {
+    if (!func) {
+      return R_NilValue;
+    }
+    return Rcpp::as<SEXP>(*this->func.get());
+  }
 
 private:
-  std::optional<Rcpp::Function> func;
+  std::shared_ptr<Rcpp::Function> func;
 };
 
-#endif // RPOET_H_
+} // namespace poet
\ No newline at end of file
diff --git a/src/Base/argh.hpp b/src/Base/argh.hpp
deleted file mode 100644
index 6b5f1f9c0..000000000
--- a/src/Base/argh.hpp
+++ /dev/null
@@ -1,459 +0,0 @@
-/*
-** Copyright (c) 2016, Adi Shavit All rights reserved.
-**
-** Redistribution and use in source and binary forms, with or without
-** modification, are permitted provided that the following conditions are met:
-**
-** * Redistributions of source code must retain the above copyright notice, this
-**   list of conditions and the following disclaimer. * Redistributions in
-**   binary form must reproduce the above copyright notice, this list of
-**   conditions and the following disclaimer in the documentation and/or other
-**   materials provided with the distribution. * Neither the name of nor the
-**   names of its contributors may be used to endorse or promote products
-**   derived from this software without specific prior written permission.
-**
-** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-** AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-** ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-** LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-** INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-** CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-** ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-** POSSIBILITY OF SUCH DAMAGE.
- */
-
-#pragma once
-
-#include <algorithm>
-#include <sstream>
-#include <string>
-#include <vector>
-#include <set>
-#include <map>
-#include <cassert>
-
-namespace argh
-{
-   // Terminology:
-   // A command line is composed of 2 types of args:
-   // 1. Positional args, i.e. free standing values
-   // 2. Options: args beginning with '-'. We identify two kinds:
-   //    2.1: Flags: boolean options =>  (exist ? true : false)
-   //    2.2: Parameters: a name followed by a non-option value
-
-#if !defined(__GNUC__) || (__GNUC__ >= 5)
-   using string_stream = std::istringstream;
-#else
-    // Until GCC 5, istringstream did not have a move constructor.
-    // stringstream_proxy is used instead, as a workaround.
-   class stringstream_proxy
-   {
-   public:
-      stringstream_proxy() = default;
-
-      // Construct with a value.
-      stringstream_proxy(std::string const& value) :
-         stream_(value)
-      {}
-
-      // Copy constructor.
-      stringstream_proxy(const stringstream_proxy& other) :
-         stream_(other.stream_.str())
-      {
-         stream_.setstate(other.stream_.rdstate());
-      }
-
-      void setstate(std::ios_base::iostate state) { stream_.setstate(state); }
-
-      // Stream out the value of the parameter.
-      // If the conversion was not possible, the stream will enter the fail state,
-      // and operator bool will return false.
-      template<typename T>
-      stringstream_proxy& operator >> (T& thing)
-      {
-         stream_ >> thing;
-         return *this;
-      }
-
-
-      // Get the string value.
-      std::string str() const { return stream_.str(); }
-
-      std::stringbuf* rdbuf() const { return stream_.rdbuf(); }
-
-      // Check the state of the stream. 
-      // False when the most recent stream operation failed
-      operator bool() const { return !!stream_; }
-
-      ~stringstream_proxy() = default;
-   private:
-      std::istringstream stream_;
-   };
-   using string_stream = stringstream_proxy;
-#endif
-
-   class parser
-   {
-   public:
-      enum Mode { PREFER_FLAG_FOR_UNREG_OPTION = 1 << 0,
-                  PREFER_PARAM_FOR_UNREG_OPTION = 1 << 1,
-                  NO_SPLIT_ON_EQUALSIGN = 1 << 2,
-                  SINGLE_DASH_IS_MULTIFLAG = 1 << 3,
-                };
-
-      parser() = default;
-
-      parser(std::initializer_list<char const* const> pre_reg_names)
-      {  add_params(pre_reg_names); }
-
-      parser(const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION)
-      {  parse(argv, mode); }
-
-      parser(int argc, const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION)
-      {  parse(argc, argv, mode); }
-
-      void add_param(std::string const& name);
-      void add_params(std::initializer_list<char const* const> init_list);
-
-      void parse(const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION);
-      void parse(int argc, const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION);
-
-      std::multiset<std::string>          const& flags()    const { return flags_;    }
-      std::map<std::string, std::string>  const& params()   const { return params_;   }
-      std::vector<std::string>            const& pos_args() const { return pos_args_; }
-
-      // begin() and end() for using range-for over positional args.
-      std::vector<std::string>::const_iterator begin() const { return pos_args_.cbegin(); }
-      std::vector<std::string>::const_iterator end()   const { return pos_args_.cend();   }
-      size_t size()                                    const { return pos_args_.size();   }
-
-      //////////////////////////////////////////////////////////////////////////
-      // Accessors
-
-      // flag (boolean) accessors: return true if the flag appeared, otherwise false.
-      bool operator[](std::string const& name) const;
-
-      // multiple flag (boolean) accessors: return true if at least one of the flag appeared, otherwise false.
-      bool operator[](std::initializer_list<char const* const> init_list) const;
-
-      // returns positional arg string by order. Like argv[] but without the options
-      std::string const& operator[](size_t ind) const;
-
-      // returns a std::istream that can be used to convert a positional arg to a typed value.
-      string_stream operator()(size_t ind) const;
-
-      // same as above, but with a default value in case the arg is missing (index out of range).
-      template<typename T>
-      string_stream operator()(size_t ind, T&& def_val) const;
-
-      // parameter accessors, give a name get an std::istream that can be used to convert to a typed value.
-      // call .str() on result to get as string
-      string_stream operator()(std::string const& name) const;
-
-      // accessor for a parameter with multiple names, give a list of names, get an std::istream that can be used to convert to a typed value.
-      // call .str() on result to get as string
-      // returns the first value in the list to be found.
-      string_stream operator()(std::initializer_list<char const* const> init_list) const;
-
-      // same as above, but with a default value in case the param was missing.
-      // Non-string def_val types must have an operator<<() (output stream operator)
-      // If T only has an input stream operator, pass the string version of the type as in "3" instead of 3.
-      template<typename T>
-      string_stream operator()(std::string const& name, T&& def_val) const;
-
-      // same as above but for a list of names. returns the first value to be found.
-      template<typename T>
-      string_stream operator()(std::initializer_list<char const* const> init_list, T&& def_val) const;
-
-   private:
-      string_stream bad_stream() const;
-      std::string trim_leading_dashes(std::string const& name) const;
-      bool is_number(std::string const& arg) const;
-      bool is_option(std::string const& arg) const;
-      bool got_flag(std::string const& name) const;
-      bool is_param(std::string const& name) const;
-
-   private:
-      std::vector<std::string> args_;
-      std::map<std::string, std::string> params_;
-      std::vector<std::string> pos_args_;
-      std::multiset<std::string> flags_;
-      std::set<std::string> registeredParams_;
-      std::string empty_;
-   };
-
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline void parser::parse(const char * const argv[], int mode)
-   {
-      int argc = 0;
-      for (auto argvp = argv; *argvp; ++argc, ++argvp);
-      parse(argc, argv, mode);
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline void parser::parse(int argc, const char* const argv[], int mode /*= PREFER_FLAG_FOR_UNREG_OPTION*/)
-   {
-      // convert to strings
-      args_.resize(argc);
-      std::transform(argv, argv + argc, args_.begin(), [](const char* const arg) { return arg;  });
-
-      // parse line
-      for (auto i = 0u; i < args_.size(); ++i)
-      {
-         if (!is_option(args_[i]))
-         {
-            pos_args_.emplace_back(args_[i]);
-            continue;
-         }
-
-         auto name = trim_leading_dashes(args_[i]);
-
-         if (!(mode & NO_SPLIT_ON_EQUALSIGN))
-         {
-            auto equalPos = name.find('=');
-            if (equalPos != std::string::npos)
-            {
-               params_.insert({ name.substr(0, equalPos), name.substr(equalPos + 1) });
-               continue;
-            }
-         }
-
-         // if the option is unregistered and should be a multi-flag
-         if (1 == (args_[i].size() - name.size()) &&         // single dash
-            argh::parser::SINGLE_DASH_IS_MULTIFLAG & mode && // multi-flag mode
-            !is_param(name))                                  // unregistered
-         {
-            std::string keep_param; 
-            
-            if (!name.empty() && is_param(std::string(1ul, name.back()))) // last char is param
-            {
-               keep_param += name.back();
-               name.resize(name.size() - 1);
-            }
-
-            for (auto const& c : name)
-            {
-               flags_.emplace(std::string{ c });
-            }
-
-            if (!keep_param.empty())
-            {
-               name = keep_param;
-            }
-            else
-            {
-               continue; // do not consider other options for this arg
-            }
-         }
-
-         // any potential option will get as its value the next arg, unless that arg is an option too
-         // in that case it will be determined a flag.
-         if (i == args_.size() - 1 || is_option(args_[i + 1]))
-         {
-            flags_.emplace(name);
-            continue;
-         }
-
-         // if 'name' is a pre-registered option, then the next arg cannot be a free parameter to it is skipped
-         // otherwise we have 2 modes:
-         // PREFER_FLAG_FOR_UNREG_OPTION: a non-registered 'name' is determined a flag. 
-         //                               The following value (the next arg) will be a free parameter.
-         //
-         // PREFER_PARAM_FOR_UNREG_OPTION: a non-registered 'name' is determined a parameter, the next arg
-         //                                will be the value of that option.
-
-         assert(!(mode & argh::parser::PREFER_FLAG_FOR_UNREG_OPTION)
-             || !(mode & argh::parser::PREFER_PARAM_FOR_UNREG_OPTION));
-
-         bool preferParam = mode & argh::parser::PREFER_PARAM_FOR_UNREG_OPTION;
-
-         if (is_param(name) || preferParam)
-         {
-            params_.insert({ name, args_[i + 1] });
-            ++i; // skip next value, it is not a free parameter
-            continue;
-         }
-         else
-         {
-            flags_.emplace(name);
-         }
-      };
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline string_stream parser::bad_stream() const
-   {
-      string_stream bad;
-      bad.setstate(std::ios_base::failbit);
-      return bad;
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool parser::is_number(std::string const& arg) const
-   {
-      // inefficient but simple way to determine if a string is a number (which can start with a '-')
-      std::istringstream istr(arg);
-      double number;
-      istr >> number;
-      return !(istr.fail() || istr.bad());
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool parser::is_option(std::string const& arg) const
-   {
-      assert(0 != arg.size());
-      if (is_number(arg))
-         return false;
-      return '-' == arg[0];
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline std::string parser::trim_leading_dashes(std::string const& name) const
-   {
-      auto pos = name.find_first_not_of('-');
-      return std::string::npos != pos ? name.substr(pos) : name;
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool argh::parser::got_flag(std::string const& name) const
-   {
-      return flags_.end() != flags_.find(trim_leading_dashes(name));
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool argh::parser::is_param(std::string const& name) const
-   {
-      return registeredParams_.count(name);
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool parser::operator[](std::string const& name) const
-   {
-      return got_flag(name);
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline bool parser::operator[](std::initializer_list<char const* const> init_list) const
-   {
-      return std::any_of(init_list.begin(), init_list.end(), [&](char const* const name) { return got_flag(name); });
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline std::string const& parser::operator[](size_t ind) const
-   {
-      if (ind < pos_args_.size())
-         return pos_args_[ind];
-      return empty_;
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline string_stream parser::operator()(std::string const& name) const
-   {
-      auto optIt = params_.find(trim_leading_dashes(name));
-      if (params_.end() != optIt)
-         return string_stream(optIt->second);
-      return bad_stream();
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline string_stream parser::operator()(std::initializer_list<char const* const> init_list) const
-   {
-      for (auto& name : init_list)
-      {
-         auto optIt = params_.find(trim_leading_dashes(name));
-         if (params_.end() != optIt)
-            return string_stream(optIt->second);
-      }
-      return bad_stream();
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   template<typename T>
-   string_stream parser::operator()(std::string const& name, T&& def_val) const
-   {
-      auto optIt = params_.find(trim_leading_dashes(name));
-      if (params_.end() != optIt)
-         return string_stream(optIt->second);
-
-      std::ostringstream ostr;
-      ostr << def_val;
-      return string_stream(ostr.str()); // use default
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   // same as above but for a list of names. returns the first value to be found.
-   template<typename T>
-   string_stream parser::operator()(std::initializer_list<char const* const> init_list, T&& def_val) const
-   {
-      for (auto& name : init_list)
-      {
-         auto optIt = params_.find(trim_leading_dashes(name));
-         if (params_.end() != optIt)
-            return string_stream(optIt->second);
-      }      
-      std::ostringstream ostr;
-      ostr << def_val;
-      return string_stream(ostr.str()); // use default
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline string_stream parser::operator()(size_t ind) const
-   {
-      if (pos_args_.size() <= ind)
-         return bad_stream();
-
-      return string_stream(pos_args_[ind]);
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   template<typename T>
-   string_stream parser::operator()(size_t ind, T&& def_val) const
-   {
-      if (pos_args_.size() <= ind)
-      {
-         std::ostringstream ostr;
-         ostr << def_val;
-         return string_stream(ostr.str());
-      }
-
-      return string_stream(pos_args_[ind]);
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline void parser::add_param(std::string const& name)
-   {
-      registeredParams_.insert(trim_leading_dashes(name));
-   }
-
-   //////////////////////////////////////////////////////////////////////////
-
-   inline void parser::add_params(std::initializer_list<char const* const> init_list)
-   {
-      for (auto& name : init_list)
-         registeredParams_.insert(trim_leading_dashes(name));
-   }
-}
-
-
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 24c10c834..4b07f355a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,6 +27,16 @@ target_link_libraries(
   PUBLIC MPI::MPI_C
 )
 
+include(FetchContent)
+FetchContent_Declare(
+  cli11
+  QUIET 
+  GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
+  GIT_TAG v2.4.2
+)
+
+FetchContent_MakeAvailable(cli11)
+
 # add_library(poetlib
 #   Base/Grid.cpp
 #   Base/SimParams.cpp
@@ -75,11 +85,11 @@ file(READ "${PROJECT_SOURCE_DIR}/R_lib/ai_surrogate_model.R" R_AI_SURROGATE_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_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime CLI11::CLI11)
 target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
 
 add_executable(poet_init initializer.cpp)
-target_link_libraries(poet_init PRIVATE POETLib RRuntime)
+target_link_libraries(poet_init PRIVATE POETLib RRuntime CLI11::CLI11)
 target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
 
 install(TARGETS poet poet_init DESTINATION bin)
diff --git a/src/Chemistry/SurrogateModels/DHT.c b/src/Chemistry/SurrogateModels/DHT.c
index b0d5b3caa..e17ff6fa0 100644
--- a/src/Chemistry/SurrogateModels/DHT.c
+++ b/src/Chemistry/SurrogateModels/DHT.c
@@ -1,4 +1,3 @@
-/// Time-stamp: "Last modified 2023-08-10 11:50:46 mluebke"
 /*
 ** Copyright (C) 2017-2021 Max Luebke (University of Potsdam)
 **
diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
index eac73b14a..095f83cfa 100644
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
+++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-11-01 10:54:45 mluebke"
 
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@@ -25,6 +24,7 @@
 #include "Init/InitialList.hpp"
 #include "Rounding.hpp"
 
+#include <Rcpp/proxy/ProtectedProxy.h>
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -267,7 +267,8 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
 
   NamedVector<double> input_nv(this->output_names, cell);
 
-  const std::vector<double> eval_vec = hooks.dht_fuzz(input_nv);
+  const std::vector<double> eval_vec =
+      Rcpp::as<std::vector<double>>(hooks.dht_fuzz(input_nv));
   assert(eval_vec.size() == this->key_count);
   LookupKey vecFuzz(this->key_count + 1, {.0});
 
diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
index fed271923..e0b1c022d 100644
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
+++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-09-08 14:43:02 mluebke"
 
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
diff --git a/src/Chemistry/SurrogateModels/HashFunctions.cpp b/src/Chemistry/SurrogateModels/HashFunctions.cpp
index d93702990..623ba9d59 100644
--- a/src/Chemistry/SurrogateModels/HashFunctions.cpp
+++ b/src/Chemistry/SurrogateModels/HashFunctions.cpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-04-24 23:20:55 mluebke"
 /*
 **-----------------------------------------------------------------------------
 ** MurmurHash2 was written by Austin Appleby, and is placed in the public
diff --git a/src/Chemistry/SurrogateModels/HashFunctions.hpp b/src/Chemistry/SurrogateModels/HashFunctions.hpp
index 49e86559a..4d9881376 100644
--- a/src/Chemistry/SurrogateModels/HashFunctions.hpp
+++ b/src/Chemistry/SurrogateModels/HashFunctions.hpp
@@ -1,4 +1,3 @@
-// //  Time-stamp: "Last modified 2023-03-31 14:59:49 mluebke"
 
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
diff --git a/src/Chemistry/SurrogateModels/Interpolation.hpp b/src/Chemistry/SurrogateModels/Interpolation.hpp
index a1886965d..b6b9a69f5 100644
--- a/src/Chemistry/SurrogateModels/Interpolation.hpp
+++ b/src/Chemistry/SurrogateModels/Interpolation.hpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-16 16:49:31 mluebke"
 
 #ifndef INTERPOLATION_H_
 #define INTERPOLATION_H_
diff --git a/src/Chemistry/SurrogateModels/InterpolationModule.cpp b/src/Chemistry/SurrogateModels/InterpolationModule.cpp
index 455c96729..a2871bb8a 100644
--- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp
+++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
 #include "Init/InitialList.hpp"
 #include "Interpolation.hpp"
 
@@ -9,6 +8,7 @@
 #include "Rounding.hpp"
 
 #include <Rcpp.h>
+#include <Rcpp/proxy/ProtectedProxy.h>
 #include <Rinternals.h>
 
 #include <algorithm>
@@ -94,7 +94,8 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
     if (hooks.interp_pre.isValid()) {
       NamedVector<double> nv_in(this->out_names, work_package.input[wp_i]);
 
-      auto rm_indices = hooks.interp_pre(nv_in, pht_result.in_values);
+      std::vector<int> rm_indices = Rcpp::as<std::vector<int>>(
+          hooks.interp_pre(nv_in, pht_result.in_values));
 
       pht_result.size -= rm_indices.size();
 
diff --git a/src/Chemistry/SurrogateModels/LookupKey.hpp b/src/Chemistry/SurrogateModels/LookupKey.hpp
index f64b92721..e6dc7e697 100644
--- a/src/Chemistry/SurrogateModels/LookupKey.hpp
+++ b/src/Chemistry/SurrogateModels/LookupKey.hpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-11 10:12:52 mluebke"
 
 #ifndef LOOKUPKEY_H_
 #define LOOKUPKEY_H_
diff --git a/src/Chemistry/SurrogateModels/ProximityHashTable.cpp b/src/Chemistry/SurrogateModels/ProximityHashTable.cpp
index f1589fa21..d63ccb7ae 100644
--- a/src/Chemistry/SurrogateModels/ProximityHashTable.cpp
+++ b/src/Chemistry/SurrogateModels/ProximityHashTable.cpp
@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-15 14:50:59 mluebke"
 #include "Interpolation.hpp"
 
 #include "DHT_Wrapper.hpp"
diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp
index 3eb2fe2fd..242929b47 100644
--- a/src/Init/DiffusionInit.cpp
+++ b/src/Init/DiffusionInit.cpp
@@ -27,9 +27,9 @@ namespace poet {
 enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC };
 
 const std::map<std::uint8_t, std::string> tug_side_mapping = {
-    {tug::BC_SIDE_RIGHT, "E"},
-    {tug::BC_SIDE_LEFT, "W"},
-    {tug::BC_SIDE_TOP, "N"},
+    {tug::BC_SIDE_RIGHT	, "E"},
+    {tug::BC_SIDE_LEFT	, "W"},
+    {tug::BC_SIDE_TOP	, "N"},
     {tug::BC_SIDE_BOTTOM, "S"}};
 
 static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
@@ -146,7 +146,7 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
         extend_initial_grid(this->initial_grid, this->transport_names);
   }
 
-  for (const auto &species : this->transport_names) {
+  for (const auto &species_name : this->transport_names) {
     Rcpp::List spec_list;
 
     for (const auto &side : tug_side_mapping) {
@@ -163,8 +163,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
       if (boundaries_list.containsElementNamed(side.second.c_str())) {
         const Rcpp::List mapping = boundaries_list[side.second];
 
-        const Rcpp::NumericVector cells = mapping["cell"];
-        const Rcpp::NumericVector values = mapping["sol_id"];
+        const Rcpp::IntegerVector cells  = mapping["cell"];   // MDL 2024-06-13 
+        const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL 
         const Rcpp::CharacterVector type_str = mapping["type"];
 
         if (cells.size() != values.size()) {
@@ -180,7 +180,7 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
           } else if (type_str[i] == "constant") {
             c_type[c_id] = tug::BC_TYPE_CONSTANT;
             c_value[c_id] = Rcpp::as<TugType>(
-                resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
+                resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i]));
           } else {
             throw std::runtime_error("Unknown boundary type");
           }
@@ -191,7 +191,7 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
           Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
     }
 
-    bound_list[species] = spec_list;
+    bound_list[species_name] = spec_list;
 
     if (inner_boundaries.size() > 0) {
       const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"];
@@ -212,10 +212,10 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
         rows.push_back(inner_row_vec[i] - 1);
         cols.push_back(inner_col_vec[i] - 1);
         c_value.push_back(Rcpp::as<TugType>(resolve_R(
-            this->phreeqc_mat, Rcpp::wrap(species), inner_pqc_id_vec[i])));
+            this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i])));
       }
 
-      inner_bound[species] =
+      inner_bound[species_name] =
           Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows),
                              Rcpp::Named("col") = Rcpp::wrap(cols),
                              Rcpp::Named("value") = Rcpp::wrap(c_value));
@@ -403,4 +403,4 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
   return diff_init;
 }
 
-} // namespace poet
\ No newline at end of file
+} // namespace poet
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index 3e6ae7654..3a3c5ea23 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -215,10 +215,10 @@ private:
 
 public:
   struct ChemistryHookFunctions {
-    RHookFunction<bool> dht_fill;
-    RHookFunction<std::vector<double>> dht_fuzz;
-    RHookFunction<std::vector<std::size_t>> interp_pre;
-    RHookFunction<bool> interp_post;
+    poet::DEFunc dht_fill;
+    poet::DEFunc dht_fuzz;
+    poet::DEFunc interp_pre;
+    poet::DEFunc interp_post;
   };
 
   struct ChemistryInit {
diff --git a/src/initializer.cpp b/src/initializer.cpp
index d9a8afadf..d2e663931 100644
--- a/src/initializer.cpp
+++ b/src/initializer.cpp
@@ -1,7 +1,7 @@
 #include "Init/InitialList.hpp"
 #include "poet.hpp"
 
-#include "Base/argh.hpp"
+#include <CLI/CLI.hpp>
 
 #include <cstdlib>
 
@@ -11,31 +11,39 @@
 #include <string>
 
 int main(int argc, char **argv) {
-
-  argh::parser cmdl({"-o", "--output"});
-
-  cmdl.parse(argc, argv);
-
-  if (cmdl[{"-h", "--help"}] || cmdl.pos_args().size() != 2) {
-    std::cout << "Usage: " << argv[0]
-              << " [-o, --output output_file]"
-                 " [-s, --setwd] "
-                 " <script.R>"
-              << std::endl;
-    return EXIT_SUCCESS;
-  }
-
+  // initialize RIinside
   RInside R(argc, argv);
 
   R.parseEvalQ(init_r_library);
+  R.parseEvalQ(kin_r_library);
 
-  std::string input_script = cmdl.pos_args()[1];
+  // parse command line arguments
+  CLI::App app{"POET Initializer - Poet R script to POET qs/rds converter"};
+
+  std::string input_script;
+  app.add_option("PoetScript.R", input_script, "POET R script to convert")
+      ->required();
+
+  std::string output_file;
+  app.add_option("-n, --name", output_file,
+                 "Name of the output file without file extension");
+
+  bool setwd;
+  app.add_flag("-s, --setwd", setwd,
+               "Set working directory to the directory of the directory of the "
+               "input script")
+      ->default_val(false);
+
+  bool asRDS;
+  app.add_flag("-r, --rds", asRDS, "Save output as .rds file instead of .qs")
+      ->default_val(false);
+
+  CLI11_PARSE(app, argc, argv);
+
+  // source the input script
   std::string normalized_path_script;
   std::string in_file_name;
 
-  std::string curr_path =
-      Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(".")));
-
   try {
     normalized_path_script =
         Rcpp::as<std::string>(Rcpp::Function("normalizePath")(input_script));
@@ -51,14 +59,20 @@ int main(int argc, char **argv) {
     return EXIT_FAILURE;
   }
 
-  std::string output_file;
+  // if output file is not specified, use the input file name
+  if (output_file.empty()) {
+    std::string curr_path =
+        Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(".")));
 
-  cmdl({"-o", "--output"},
-       curr_path + "/" +
-           in_file_name.substr(0, in_file_name.find_last_of('.')) + ".rds") >>
-      output_file;
+    output_file = curr_path + "/" +
+                  in_file_name.substr(0, in_file_name.find_last_of('.'));
+  }
 
-  if (cmdl[{"-s", "--setwd"}]) {
+  // append the correct file extension
+  output_file += asRDS ? ".rds" : ".qs";
+
+  // set working directory to the directory of the input script
+  if (setwd) {
     const std::string dir_path = Rcpp::as<std::string>(
         Rcpp::Function("dirname")(normalized_path_script));
 
@@ -71,13 +85,12 @@ int main(int argc, char **argv) {
 
   init.initializeFromList(setup);
 
-  // replace file extension by .rds
-  Rcpp::Function save("saveRDS");
-
-  save(init.exportList(), Rcpp::wrap(output_file));
+  // use the generic handler defined in kin_r_library.R
+  Rcpp::Function SaveRObj_R("SaveRObj");
+  SaveRObj_R(init.exportList(), Rcpp::wrap(output_file));
 
   std::cout << "Saved result to " << output_file << std::endl;
 
   //   parseGrid(R, grid, results);
   return EXIT_SUCCESS;
-}
\ No newline at end of file
+}
diff --git a/src/poet.cpp b/src/poet.cpp
index 4a0abc2c1..9fbf94c18 100644
--- a/src/poet.cpp
+++ b/src/poet.cpp
@@ -4,7 +4,8 @@
 **
 ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
 **
-** Copyright (C) 2023-2024 Max Luebke (University of 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
@@ -22,6 +23,7 @@
 
 #include "Base/Macros.hpp"
 #include "Base/RInsidePOET.hpp"
+#include "CLI/CLI.hpp"
 #include "Chemistry/ChemistryModule.hpp"
 #include "DataStructures/Field.hpp"
 #include "Init/InitialList.hpp"
@@ -36,10 +38,9 @@
 #include <cstdlib>
 #include <memory>
 #include <mpi.h>
-#include <optional>
 #include <string>
 
-#include "Base/argh.hpp"
+#include <CLI/CLI.hpp>
 
 #include <poet.hpp>
 #include <vector>
@@ -52,17 +53,23 @@ static int MY_RANK = 0;
 
 static std::unique_ptr<Rcpp::List> global_rt_setup;
 
-// we need some layz evaluation, as we can't define the functions before the R
-// runtime is initialized
-static std::optional<Rcpp::Function> master_init_R;
-static std::optional<Rcpp::Function> master_iteration_end_R;
-static std::optional<Rcpp::Function> store_setup_R;
+// we need some lazy evaluation, as we can't define the functions
+// before the R runtime is initialized
+static poet::DEFunc master_init_R;
+static poet::DEFunc master_iteration_end_R;
+static poet::DEFunc store_setup_R;
+static poet::DEFunc ReadRObj_R;
+static poet::DEFunc SaveRObj_R;
+static poet::DEFunc source_R;
 
 static void init_global_functions(RInside &R) {
   R.parseEval(kin_r_library);
-  master_init_R = Rcpp::Function("master_init");
-  master_iteration_end_R = Rcpp::Function("master_iteration_end");
-  store_setup_R = Rcpp::Function("StoreSetup");
+  master_init_R = DEFunc("master_init");
+  master_iteration_end_R = DEFunc("master_iteration_end");
+  store_setup_R = DEFunc("StoreSetup");
+  source_R = DEFunc("source");
+  ReadRObj_R = DEFunc("ReadRObj");
+  SaveRObj_R = DEFunc("SaveRObj");
 }
 
 // HACK: this is a step back as the order and also the count of fields is
@@ -80,78 +87,98 @@ static void init_global_functions(RInside &R) {
 
 enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
 
-ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
-  // initialize argh object
-  argh::parser cmdl(argv);
+int parseInitValues(int argc, char **argv, RuntimeParameters &params) {
 
-  // if user asked for help
-  if (cmdl[{"help", "h"}]) {
-    if (MY_RANK == 0) {
-      MSG("Todo");
-      MSG("See README.md for further information.");
-    }
+  CLI::App app{"POET - Potsdam rEactive Transport simulator"};
 
-    return ParseRet::PARSER_HELP;
-  }
-  // if positional arguments are missing
-  if (!cmdl(3)) {
-    if (MY_RANK == 0) {
-      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;
-  }
-
-  // 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;
-    }
-  }
-
-  // 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;
-    }
-  }
-
-  params.print_progressbar = cmdl[{"P", "progress"}];
+  app.add_flag("-P,--progress", params.print_progress,
+               "Print progress bar during chemical simulation");
 
   /*Parse work package size*/
-  cmdl("work-package-size", CHEM_DEFAULT_WORK_PACKAGE_SIZE) >>
-      params.work_package_size;
+  app.add_option(
+         "-w,--work-package-size", params.work_package_size,
+         "Work package size to distribute to each worker for chemistry module")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::WORK_PACKAGE_SIZE_DEFAULT);
 
   /* Parse DHT arguments */
-  params.use_dht = cmdl["dht"];
-  params.use_interp = cmdl["interp"];
+  auto *dht_group = app.add_option_group("DHT", "DHT related options");
+
+  dht_group->add_flag("--dht", params.use_dht, "Enable DHT");
+
   // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
 
-  cmdl("dht-size", CHEM_DHT_SIZE_PER_PROCESS_MB) >> params.dht_size;
+  dht_group
+      ->add_option("--dht-size", params.dht_size,
+                   "DHT size per process in Megabyte")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::DHT_SIZE_DEFAULT);
   // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
   // endl;
 
-  cmdl("dht-snaps", 0) >> params.dht_snaps;
+  dht_group->add_option(
+      "--dht-snaps", params.dht_snaps,
+      "Save snapshots of DHT to disk: \n0 = disabled (default)\n1 = After "
+      "simulation\n2 = After each iteration");
 
-  params.use_interp = cmdl["interp"];
-  cmdl("interp-size", 100) >> params.interp_size;
-  cmdl("interp-min", 5) >> params.interp_min_entries;
-  cmdl("interp-bucket-entries", 20) >> params.interp_bucket_entries;
+  auto *interp_group =
+      app.add_option_group("Interpolation", "Interpolation related options");
 
-  params.use_ai_surrogate = cmdl["ai-surrogate"];
+  interp_group->add_flag("--interp", params.use_interp, "Enable interpolation");
+  interp_group
+      ->add_option("--interp-size", params.interp_size,
+                   "Size of the interpolation table in Megabyte")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_SIZE_DEFAULT);
+  interp_group
+      ->add_option("--interp-min", params.interp_min_entries,
+                   "Minimum number of entries in the interpolation table")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_MIN_ENTRIES_DEFAULT);
+  interp_group
+      ->add_option(
+          "--interp-bucket-entries", params.interp_bucket_entries,
+          "Maximum number of entries in each bucket of the interpolation table")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_BUCKET_ENTRIES_DEFAULT);
+
+  app.add_flag("--ai-surrogate", params.use_ai_surrogate,
+               "Enable AI surrogate for chemistry module");
+
+  app.add_flag("--rds", params.as_rds,
+               "Save output as .rds file instead of .qs");
+
+  std::string init_file;
+  std::string runtime_file;
+
+  app.add_option("runtime_file", runtime_file,
+                 "Runtime R script defining the simulation")
+      ->required()
+      ->check(CLI::ExistingFile);
+
+  app.add_option(
+         "init_file", init_file,
+         "Initial R script defining the simulation, produced by poet_init")
+      ->required()
+      ->check(CLI::ExistingFile);
+
+  app.add_option("out_dir", params.out_dir,
+                 "Output directory of the simulation")
+      ->required();
+
+  try {
+    app.parse(argc, argv);
+  } catch (const CLI::ParseError &e) {
+    app.exit(e);
+    return -1;
+  }
+
+  // set the output extension
+  params.out_ext = params.as_rds ? "rds" : "qs";
 
   if (MY_RANK == 0) {
     // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
+    MSG("Output format/extension is " + params.out_ext);
     MSG("Work Package Size: " + std::to_string(params.work_package_size));
     MSG("DHT is " + BOOL_PRINT(params.use_dht));
     MSG("AI Surrogate is " + BOOL_PRINT(params.use_ai_surrogate));
@@ -178,14 +205,6 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
           std::to_string(params.interp_bucket_entries));
     }
   }
-
-  std::string init_file;
-  std::string runtime_file;
-
-  cmdl(1) >> runtime_file;
-  cmdl(2) >> init_file;
-  cmdl(3) >> params.out_dir;
-
   // chem_params.dht_outdir = out_dir;
 
   /* distribute information to R runtime */
@@ -207,16 +226,17 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
   // R["dht_log"] = simparams.dht_log;
 
   try {
-    Rcpp::Function source("source");
-    Rcpp::Function readRDS("readRDS");
 
-    Rcpp::List init_params_ = readRDS(init_file);
+    Rcpp::List init_params_(ReadRObj_R(init_file));
     params.init_params = init_params_;
 
     global_rt_setup = std::make_unique<Rcpp::List>();
-    *global_rt_setup = source(runtime_file, Rcpp::Named("local", true));
+    *global_rt_setup = source_R(runtime_file, Rcpp::Named("local", true));
     *global_rt_setup = global_rt_setup->operator[]("value");
 
+    // MDL add "out_ext" for output format to R setup
+    (*global_rt_setup)["out_ext"] = params.out_ext;
+
     params.timesteps =
         Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("timesteps"));
 
@@ -255,11 +275,11 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
    * master for the following loop) */
   uint32_t maxiter = params.timesteps.size();
 
-  if (params.print_progressbar) {
+  if (params.print_progress) {
     chem.setProgressBarPrintout(true);
   }
   R["TMP_PROPS"] = Rcpp::wrap(chem.getField().GetProps());
-  
+
   /* SIMULATION LOOP */
   double dSimTime{0};
   for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
@@ -283,18 +303,20 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
       double ai_start_t = MPI_Wtime();
       // Save current values from the tug field as predictor for the ai step
       R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
-      R.parseEval(std::string(
-        "predictors <- setNames(data.frame(matrix(TMP, nrow=" +
-        std::to_string(chem.getField().GetRequestedVecSize()) + ")), TMP_PROPS)"));
+      R.parseEval(
+          std::string("predictors <- setNames(data.frame(matrix(TMP, nrow=" +
+                      std::to_string(chem.getField().GetRequestedVecSize()) +
+                      ")), TMP_PROPS)"));
       R.parseEval("predictors <- predictors[ai_surrogate_species]");
 
-      // Apply preprocessing 
+      // Apply preprocessing
       MSG("AI Preprocessing");
       R.parseEval("predictors_scaled <- preprocess(predictors)");
 
       // Predict
       MSG("AI Predict");
-      R.parseEval("aipreds_scaled <- prediction_step(model, predictors_scaled)");
+      R.parseEval(
+          "aipreds_scaled <- prediction_step(model, predictors_scaled)");
 
       // Apply postprocessing
       MSG("AI Postprocesing");
@@ -302,24 +324,26 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
 
       // Validate prediction and write valid predictions to chem field
       MSG("AI Validate");
-      R.parseEval("validity_vector <- validate_predictions(predictors, aipreds)");
-      
+      R.parseEval(
+          "validity_vector <- validate_predictions(predictors, aipreds)");
+
       MSG("AI Marking accepted");
       chem.set_ai_surrogate_validity_vector(R.parseEval("validity_vector"));
 
       MSG("AI TempField");
-      std::vector<std::vector<double>> RTempField = R.parseEval("set_valid_predictions(predictors,\
+      std::vector<std::vector<double>> RTempField =
+          R.parseEval("set_valid_predictions(predictors,\
                                                                                        aipreds,\
                                                                                        validity_vector)");
 
       MSG("AI Set Field");
-      Field predictions_field = Field(R.parseEval("nrow(predictors)"), 
-                                      RTempField,
-                                      R.parseEval("colnames(predictors)"));
+      Field predictions_field =
+          Field(R.parseEval("nrow(predictors)"), RTempField,
+                R.parseEval("colnames(predictors)"));
 
       MSG("AI Update");
       chem.getField().update(predictions_field);
-      double ai_end_t = MPI_Wtime();  
+      double ai_end_t = MPI_Wtime();
       R["ai_prediction_time"] = ai_end_t - ai_start_t;
     }
 
@@ -330,16 +354,18 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
       double ai_start_t = MPI_Wtime();
 
       R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
-      R.parseEval(std::string(
-        "targets <- setNames(data.frame(matrix(TMP, nrow=" +
-        std::to_string(chem.getField().GetRequestedVecSize()) + ")), TMP_PROPS)"));
+      R.parseEval(
+          std::string("targets <- setNames(data.frame(matrix(TMP, nrow=" +
+                      std::to_string(chem.getField().GetRequestedVecSize()) +
+                      ")), TMP_PROPS)"));
       R.parseEval("targets <- targets[ai_surrogate_species]");
-      
+
       // TODO: Check how to get the correct columns
       R.parseEval("target_scaled <- preprocess(targets)");
 
       MSG("AI: incremental training");
-      R.parseEval("model <- training_step(model, predictors_scaled, target_scaled, validity_vector)");
+      R.parseEval("model <- training_step(model, predictors_scaled, "
+                  "target_scaled, validity_vector)");
       double ai_end_t = MPI_Wtime();
       R["ai_training_time"] = ai_end_t - ai_start_t;
     }
@@ -463,17 +489,24 @@ int main(int argc, char *argv[]) {
       MSG("Running POET version " + std::string(poet_version));
     }
 
+    init_global_functions(R);
+
     RuntimeParameters run_params;
 
-    switch (parseInitValues(argv, run_params)) {
-    case ParseRet::PARSER_ERROR:
-    case ParseRet::PARSER_HELP:
+    if (parseInitValues(argc, argv, run_params) != 0) {
       MPI_Finalize();
       return 0;
-    case ParseRet::PARSER_OK:
-      break;
     }
 
+    // switch (parseInitValues(argc, argv, run_params)) {
+    // case ParseRet::PARSER_ERROR:
+    // case ParseRet::PARSER_HELP:
+    //   MPI_Finalize();
+    //   return 0;
+    // case ParseRet::PARSER_OK:
+    //   break;
+    // }
+
     InitialList init_list(R);
     init_list.importList(run_params.init_params, MY_RANK != 0);
 
@@ -487,6 +520,7 @@ int main(int argc, char *argv[]) {
                               init_list.getChemistryInit(), MPI_COMM_WORLD);
 
     const ChemistryModule::SurrogateSetup surr_setup = {
+
         getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
         run_params.use_dht,
         run_params.dht_size,
@@ -501,33 +535,37 @@ int main(int argc, char *argv[]) {
     if (MY_RANK > 0) {
       chemistry.WorkerLoop();
     } else {
-      init_global_functions(R);
       // R.parseEvalQ("mysetup <- setup");
       // // if (MY_RANK == 0) { // get timestep vector from
       // // grid_init function ... //
-      *global_rt_setup =
-          master_init_R.value()(*global_rt_setup, run_params.out_dir,
-                                init_list.getInitialGrid().asSEXP());
+
+      *global_rt_setup = master_init_R(*global_rt_setup, run_params.out_dir,
+                                       init_list.getInitialGrid().asSEXP());
+
       // MDL: store all parameters
       // MSG("Calling R Function to store calling parameters");
       // R.parseEvalQ("StoreSetup(setup=mysetup)");
+      R["out_ext"] = run_params.out_ext;
+      R["out_dir"] = run_params.out_dir;
+
       if (run_params.use_ai_surrogate) {
         /* Incorporate ai surrogate from R */
         R.parseEvalQ(ai_surrogate_r_library);
         /* Use dht species for model input and output */
-        R["ai_surrogate_species"] = init_list.getChemistryInit().dht_species.getNames();
-        R["out_dir"] = run_params.out_dir;
-        
-        const std::string ai_surrogate_input_script = init_list.getChemistryInit().ai_surrogate_input_script;
+        R["ai_surrogate_species"] =
+            init_list.getChemistryInit().dht_species.getNames();
+
+        const std::string ai_surrogate_input_script =
+            init_list.getChemistryInit().ai_surrogate_input_script;
+
+        MSG("AI: sourcing user-provided script");
+        R.parseEvalQ(ai_surrogate_input_script);
 
-    	  MSG("AI: sourcing user-provided script");
-	      R.parseEvalQ(ai_surrogate_input_script);
-	      
         MSG("AI: initialize AI model");
-	      R.parseEval("model <- initiate_model()");
+        R.parseEval("model <- initiate_model()");
         R.parseEval("gpu_info()");
-        }
-        
+      }
+
       MSG("Init done on process with rank " + std::to_string(MY_RANK));
 
       // MPI_Barrier(MPI_COMM_WORLD);
@@ -543,14 +581,15 @@ int main(int argc, char *argv[]) {
 
       R["profiling"] = profiling;
       R["setup"] = *global_rt_setup;
+      R["setup$out_ext"] = run_params.out_ext;
 
       string r_vis_code;
-      r_vis_code =
-	"saveRDS(profiling, file=paste0(setup$out_dir,'/timings.rds'));";
+      r_vis_code = "SaveRObj(x = profiling, path = paste0(out_dir, "
+                   "'/timings.', setup$out_ext));";
       R.parseEval(r_vis_code);
 
       MSG("Done! Results are stored as R objects into <" + run_params.out_dir +
-          "/timings.rds>");
+          "/timings." + run_params.out_ext);
     }
   }
 
diff --git a/src/poet.hpp.in b/src/poet.hpp.in
index cca89e264..0e2409f87 100644
--- a/src/poet.hpp.in
+++ b/src/poet.hpp.in
@@ -20,7 +20,7 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
-#pragma once  
+#pragma once
 
 #include <cstdint>
 #include <set>
@@ -35,62 +35,37 @@ static const char *poet_version = "@POET_VERSION@";
 static const inline std::string kin_r_library = R"(@R_KIN_LIB@)";
 
 static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
-static const inline std::string ai_surrogate_r_library = R"(@R_AI_SURROGATE_LIB@)";
+static const inline std::string ai_surrogate_r_library =
+    R"(@R_AI_SURROGATE_LIB@)";
 static const inline std::string r_runtime_parameters = "mysetup";
 
-const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
-                                     "interp", "ai-surrogate"};
-const std::set<std::string> paramlist{
-    "work-package-size", "dht-strategy", "dht-size",   "dht-snaps",
-    "dht-file",          "interp-size",  "interp-min", "interp-bucket-entries"};
-
-constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32;
-
-constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3;
-
 struct RuntimeParameters {
   std::string out_dir;
   std::vector<double> timesteps;
 
-  bool print_progressbar;
-  uint32_t work_package_size;
-
   Rcpp::List init_params;
 
-  bool use_dht;
+  bool as_rds = false;
+  std::string out_ext; // MDL added to accomodate for qs::qsave/qread
+
+  bool print_progress = false;
+
+  static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
+  std::uint32_t work_package_size;
+
+  bool use_dht = false;
+  static constexpr std::uint32_t DHT_SIZE_DEFAULT = 1.5E3;
   std::uint32_t dht_size;
+  static constexpr std::uint8_t DHT_SNAPS_DEFAULT = 0;
   std::uint8_t dht_snaps;
 
-  bool use_interp;
+  bool use_interp = false;
+  static constexpr std::uint32_t INTERP_SIZE_DEFAULT = 100;
   std::uint32_t interp_size;
+  static constexpr std::uint32_t INTERP_MIN_ENTRIES_DEFAULT = 5;
   std::uint32_t interp_min_entries;
+  static constexpr std::uint32_t INTERP_BUCKET_ENTRIES_DEFAULT = 20;
   std::uint32_t interp_bucket_entries;
 
-  bool use_ai_surrogate;
-  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);
-  };
+  bool use_ai_surrogate = false;
 };
diff --git a/test/testField.cpp b/test/testField.cpp
index 0800b4dbd..51858ecc4 100644
--- a/test/testField.cpp
+++ b/test/testField.cpp
@@ -89,14 +89,14 @@ TEST_CASE("Field") {
   }
 
   SUBCASE("Apply R function (set Na to zero)") {
-    RHookFunction<Field> to_call(R, "simple_field");
+    poet::DEFunc to_call("simple_field");
     Field field_proc = to_call(dut.asSEXP());
 
     CHECK_EQ(field_proc["Na"], FieldColumn(dut.GetRequestedVecSize(), 0));
   }
 
   SUBCASE("Apply R function (add two fields)") {
-    RHookFunction<Field> to_call(R, "extended_field");
+    poet::DEFunc to_call("extended_field");
     Field field_proc = to_call(dut.asSEXP(), dut.asSEXP());
 
     CHECK_EQ(field_proc["Na"],
diff --git a/test/testNamedVector.cpp b/test/testNamedVector.cpp
index 7b86c7496..71d575ba0 100644
--- a/test/testNamedVector.cpp
+++ b/test/testNamedVector.cpp
@@ -9,7 +9,7 @@
 #include "testDataStructures.hpp"
 
 TEST_CASE("NamedVector") {
-  RInsidePOET &R = RInsidePOET::getInstance();
+  poet::RInsidePOET &R = poet::RInsidePOET::getInstance();
 
   R["sourcefile"] = RInside_source_file;
   R.parseEval("source(sourcefile)");
@@ -36,14 +36,14 @@ TEST_CASE("NamedVector") {
   }
 
   SUBCASE("Apply R function (set to zero)") {
-    RHookFunction<poet::NamedVector<double>> to_call(R, "simple_named_vec");
+    poet::DEFunc to_call("simple_named_vec");
     nv = to_call(nv);
 
     CHECK_EQ(nv[2], 0);
   }
 
   SUBCASE("Apply R function (second NamedVector)") {
-    RHookFunction<poet::NamedVector<double>> to_call(R, "extended_named_vec");
+    poet::DEFunc to_call("extended_named_vec");
 
     const std::vector<std::string> names{{"C", "H", "Mg"}};
     const std::vector<double> values{{0, 1, 2}};
@@ -56,8 +56,8 @@ TEST_CASE("NamedVector") {
   }
 
   SUBCASE("Apply R function (check if zero)") {
-    RHookFunction<bool> to_call(R, "bool_named_vec");
+    poet::DEFunc to_call("bool_named_vec");
 
-    CHECK_FALSE(to_call(nv));
+    CHECK_FALSE(Rcpp::as<bool>(to_call(nv)));
   }
 }