diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R
index f2dfbb71c..83eaa7293 100644
--- a/R_lib/init_r_lib.R
+++ b/R_lib/init_r_lib.R
@@ -34,4 +34,48 @@ pqc_to_grid <- function(pqc_in, grid) {
 pqc_init <- pqc_to_grid(input, grid_def)
 test <- pqc_init
 
+modify_module_sizes <- function(mod_sizes, pqc_mat, init_grid) {
+    # Find all unique IDs in init_grid
+    unique_ids <- unique(as.vector(init_grid$id))
+
+    # remove rows from pqc_mat that are not in init_grid
+    pqc_mat <- as.data.frame(pqc_mat)
+    pqc_mat <- pqc_mat[pqc_mat$id %in% unique_ids, ]
+
+    # Find the column indices where all rows are NA
+    na_cols <- which(colSums(is.na(pqc_mat)) == nrow(pqc_mat))
+
+    # Build cumsum over mod_sizes
+    cum_mod_sizes <- cumsum(mod_sizes)
+
+    # Find the indices where the value of na_cols is equal to the value of cum_mod_sizes
+    idx <- which(cum_mod_sizes %in% na_cols)
+
+    # Set the value of mod_sizes to 0 at the indices found in the previous step
+    mod_sizes[idx] <- 0
+
+    return(mod_sizes)
+}
+
+# mod_sizes <- c(7, 0, 4, 2, 0)
+
+# unique_ids <- unique(as.vector(pqc_init$id))
+
+# tmp <- as.data.frame(input)
+# pqc_test <- tmp[tmp$id %in% unique_ids, ]
+# na_cols <- which(colSums(is.na(pqc_test)) == nrow(pqc_test))
+# cum_mod_sizes <- cumsum(mod_sizes)
+
+# # Get the indices of the columns of cum_mod_sizes where the value of the column is equal to the value of na_cols
+# idx <- which(cum_mod_sizes %in% na_cols)
+# mod_sizes[idx] <- 0
+
+# idx <- which(na_cols == cum_mod_sizes)
+
+
+
+# idx <- which(na_cols[1] >= cum_mod_sizes)
+
+# mod_sizes <- modify_module_sizes(mod_sizes, pqc_init, pqc_init)
+
 # remove column with all NA
diff --git a/ext/iphreeqc b/ext/iphreeqc
index f3f86bb4a..54136c8e0 160000
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@@ -1 +1 @@
-Subproject commit f3f86bb4ac1d3b70aec7437781875072f2896ae3
+Subproject commit 54136c8e0c638d0dc47e85cef1428822b6c1512d
diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp
index a45c1f428..bb29615b2 100644
--- a/src/Init/GridInit.cpp
+++ b/src/Init/GridInit.cpp
@@ -3,9 +3,12 @@
 #include <IPhreeqcPOET.hpp>
 #include <Rcpp/Function.h>
 #include <Rcpp/vector/instantiation.h>
+#include <algorithm>
+#include <cstdint>
 #include <map>
 #include <regex>
 #include <string>
+#include <vector>
 
 namespace poet {
 static Rcpp::NumericMatrix pqcScriptToGrid(IPhreeqcPOET &phreeqc, RInside &R) {
@@ -55,6 +58,21 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
   return raws;
 }
 
+static inline IPhreeqcPOET::ModulesArray
+modifyModuleSizes(IPhreeqcPOET::ModulesArray sizes,
+                  const Rcpp::NumericMatrix &phreeqc_mat,
+                  const Rcpp::List &initial_grid) {
+  std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end());
+
+  Rcpp::Function modify_sizes("modify_module_sizes");
+  sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
+      modify_sizes(sizes_vec, phreeqc_mat, initial_grid));
+
+  std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin());
+
+  return sizes;
+}
+
 void InitialList::initGrid(const Rcpp::List &grid_input) {
   // parse input values
   const std::string script = Rcpp::as<std::string>(grid_input["pqc_in"]);
@@ -105,6 +123,14 @@ void InitialList::initGrid(const Rcpp::List &grid_input) {
   this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
   this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
 
+  this->module_sizes = modifyModuleSizes(phreeqc.getModuleSizes(),
+                                         this->phreeqc_mat, this->initial_grid);
+
+  // print module sizes
+  for (std::size_t i = 0; i < this->module_sizes.size(); i++) {
+    std::cout << this->module_sizes[i] << std::endl;
+  }
+
   this->pqc_raw_dumps = replaceRawKeywordIDs(phreeqc.raw_dumps());
 
   R["pqc_mat"] = this->phreeqc_mat;
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index e3f86ebf6..d236d19be 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -1,12 +1,12 @@
 #pragma once
 
-#include "../DataStructures/Field.hpp"
-
 #include <RInside.h>
 #include <Rcpp.h>
 #include <Rcpp/vector/instantiation.h>
 #include <cstdint>
 
+#include <IPhreeqcPOET.hpp>
+
 namespace poet {
 
 class InitialList {
@@ -40,5 +40,8 @@ private:
 
   // Initialized by grid, modified by chemistry
   std::map<int, std::string> pqc_raw_dumps;
+
+  // Chemistry members
+  IPhreeqcPOET::ModulesArray module_sizes;
 };
 } // namespace poet
\ No newline at end of file