diff --git a/include/poet/DiffusionModule.hpp b/include/poet/DiffusionModule.hpp
index 05b48172f..ce486b9e2 100644
--- a/include/poet/DiffusionModule.hpp
+++ b/include/poet/DiffusionModule.hpp
@@ -24,6 +24,7 @@
 #include "poet/SimParams.hpp"
 #include <array>
 #include <bits/stdint-uintn.h>
+#include <cmath>
 #include <poet/Grid.hpp>
 #include <string>
 #include <tug/BoundaryCondition.hpp>
@@ -85,6 +86,8 @@ private:
 
   void initialize(poet::DiffusionParams args);
 
+  void RoundToZero(double *field, uint32_t cell_count) const;
+
   Grid &grid;
   uint8_t dim;
 
diff --git a/src/DiffusionModule.cpp b/src/DiffusionModule.cpp
index bf8052427..f6282572e 100644
--- a/src/DiffusionModule.cpp
+++ b/src/DiffusionModule.cpp
@@ -23,6 +23,7 @@
 #include "tug/Diffusion.hpp"
 #include <Rcpp.h>
 #include <algorithm>
+#include <bits/stdint-intn.h>
 #include <cstdint>
 #include <poet/ChemSimSeq.hpp>
 #include <poet/DiffusionModule.hpp>
@@ -36,6 +37,8 @@
 
 using namespace poet;
 
+static constexpr double ZERO_MULTIPLICATOR = 10E-14;
+
 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};
@@ -180,6 +183,11 @@ void DiffusionModule::simulate(double dt) {
     } else {
       tug::diffusion::ADI_2D(this->diff_input, in_field, in_alpha.data());
     }
+
+    // TODO: do not use hardcoded index for O, H and charge
+    if (i > 2) {
+      this->RoundToZero(in_field, this->n_cells_per_prop);
+    }
   }
 
   std::cout << " done!\n";
@@ -188,6 +196,12 @@ void DiffusionModule::simulate(double dt) {
 
   transport_t += sim_a_transport - sim_b_transport;
 }
+inline void DiffusionModule::RoundToZero(double *field,
+                                         uint32_t cell_count) const {
+  for (uint32_t i = 0; i < cell_count; i++) {
+    field[i] = ((int32_t)(field[i] / ZERO_MULTIPLICATOR)) * ZERO_MULTIPLICATOR;
+  }
+}
 
 void DiffusionModule::end() {
   // R["simtime_transport"] = transport_t;