diff --git a/runtime/boundary_weights.hpp b/runtime/boundary_weights.hpp
index 59c57ca0d2090d2c8ff5f534df9d4201bca4da10..3678c899f46c1ae75820545016beb2e72ca611f6 100644
--- a/runtime/boundary_weights.hpp
+++ b/runtime/boundary_weights.hpp
@@ -51,7 +51,7 @@ void compute_boundary_weights(
         if( pos_x > xmin && pos_x <= xmax &&
             pos_y > ymin && pos_y <= ymax &&
             pos_z > zmin && pos_z <= zmax) {
-                *comp_weight++;
+                (*comp_weight)++;
         }
     }
 
@@ -63,7 +63,7 @@ void compute_boundary_weights(
         if( pos_x > xmin && pos_x <= xmax &&
             pos_y > ymin && pos_y <= ymax &&
             pos_z > zmin && pos_z <= zmax) {
-                *comm_weight++;
+                (*comm_weight)++;
         }
     }
     #else
diff --git a/runtime/domain/ParticleDataHandling.hpp b/runtime/domain/ParticleDataHandling.hpp
index 502b9a77a5efb6bbeffa5a9d8fbe36ac6055843b..490b0b39b897491e3b71cb5f31e471b9be52ea9b 100644
--- a/runtime/domain/ParticleDataHandling.hpp
+++ b/runtime/domain/ParticleDataHandling.hpp
@@ -1,6 +1,12 @@
 #include <blockforest/BlockForest.h>
 #include <blockforest/BlockDataHandling.h>
 
+namespace pairs {
+
+class PairsSimulation;
+
+}
+
 namespace walberla {
 
 namespace internal {
@@ -22,53 +28,51 @@ inline bool operator==(const ParticleDeleter& lhs, const ParticleDeleter& rhs) {
 
 } // namespace internal
 
-class ParticleDataHandling : public blockforest::BlockDataHandling<internal::ParticleDeleter>{
+class ParticleDataHandling : public blockforest::BlockDataHandling<internal::ParticleDeleter> {
 private:
     pairs::PairsSimulation *ps;
 
 public:
-    ParticleDataHandling(pairs::PairsSimulation *_ps) { ps = _ps; }
-    virtual ~ParticleDataHandling() {}
+    ParticleDataHandling(pairs::PairsSimulation *ps_) : ps(ps_) {}
+    ~ParticleDataHandling() override = default;
 
-    virtual internal::ParticleDeleter *initialize(IBlock *const block) WALBERLA_OVERRIDE {
+    internal::ParticleDeleter *initialize(IBlock *const block) override {
         return new internal::ParticleDeleter(block->getAABB());
     }
 
-    virtual void serialize(IBlock *const block, const BlockDataID& id, mpi::SendBuffer& buffer) WALBERLA_OVERRIDE {
+    void serialize(IBlock *const block, const BlockDataID& id, mpi::SendBuffer& buffer) override {
         serializeImpl(static_cast<Block *const>(block), id, buffer, 0, false);
     }
 
-    virtual internal::ParticleDeleter* deserialize(IBlock *const block) WALBERLA_OVERRIDE {
+    internal::ParticleDeleter* deserialize(IBlock *const block) override {
         return initialize(block);
     }
 
-    virtual void deserialize(IBlock *const block, const BlockDataID& id, mpi::RecvBuffer& buffer) WALBERLA_OVERRIDE {
+    void deserialize(IBlock *const block, const BlockDataID& id, mpi::RecvBuffer& buffer) override {
         deserializeImpl(block, id, buffer);
     }
 
-    virtual void serializeCoarseToFine(Block *const block, const BlockDataID& id, mpi::SendBuffer& buffer, const uint_t child)
-        WALBERLA_OVERRIDE {
+    void serializeCoarseToFine(Block *const block, const BlockDataID& id, mpi::SendBuffer& buffer, const uint_t child) override {
         serializeImpl(block, id, buffer, child, true);
     }
 
-    virtual void serializeFineToCoarse(Block *const block, const BlockDataID& id, mpi::SendBuffer& buffer) WALBERLA_OVERRIDE {
+    void serializeFineToCoarse(Block *const block, const BlockDataID& id, mpi::SendBuffer& buffer) override {
         serializeImpl(block, id, buffer, 0, false);
     }
 
-    virtual internal::ParticleDeleter *deserializeCoarseToFine(Block *const block) WALBERLA_OVERRIDE {
+    internal::ParticleDeleter *deserializeCoarseToFine(Block *const block) override {
         return initialize(block);
     }
 
-    virtual internal::ParticleDeleter *deserializeFineToCoarse(Block *const block) WALBERLA_OVERRIDE {
+    internal::ParticleDeleter *deserializeFineToCoarse(Block *const block) override {
         return initialize(block);
     }
 
-    virtual void deserializeCoarseToFine(Block *const block, const BlockDataID& id, mpi::RecvBuffer& buffer) WALBERLA_OVERRIDE {
+    void deserializeCoarseToFine(Block *const block, const BlockDataID& id, mpi::RecvBuffer& buffer) override {
         deserializeImpl(block, id, buffer);
     }
 
-    virtual void deserializeFineToCoarse(Block *const block, const BlockDataID& id, mpi::RecvBuffer& buffer, const uint_t)
-        WALBERLA_OVERRIDE {
+    void deserializeFineToCoarse(Block *const block, const BlockDataID& id, mpi::RecvBuffer& buffer, const uint_t) override {
         deserializeImpl(block, id, buffer);
     }
 
@@ -95,11 +99,13 @@ public:
             aabb_check[5] = aabb.zMax();
         }
 
-        for(auto& p: ps->getNonVolatileProperties()) {
-            ps->copyPropertyToHost(p, WriteAfterRead);
+        for(auto& prop: ps->getProperties()) {
+            if(!prop.isVolatile()) {
+                ps->copyPropertyToHost(prop, WriteAfterRead);
+            }
         }
 
-        auto position = ps->getPropertyByName("position");
+        auto position = ps->getAsVectorProperty(ps->getPropertyByName("position"));
         int nlocal = ps->getTrackedVariableAsInteger("nlocal");
         int i = 0;
         int nserialized = 0;
@@ -115,51 +121,52 @@ public:
 
                 nlocal--;
 
-                for(auto &p: ps->getNonVolatileProperties()) {
-                    auto prop = ps->getProperty(p_id);
-                    auto prop_type = prop.getType();
-
-                    if(prop_type == Prop_Vector) {
-                        auto vector_ptr = ps->getAsVectorProperty(prop);
-                        constexpr int nelems = 3;
-
-                        for(int e = 0; e < nelems; e++) {
-                            buffer << vector_ptr(i, e);
-                            vector_ptr(i, e) = vector_ptr(nlocal, e);
+                for(auto &prop: ps->getProperties()) {
+                    if(!prop.isVolatile()) {
+                        auto prop_type = prop.getType();
+
+                        if(prop_type == Prop_Vector) {
+                            auto vector_ptr = ps->getAsVectorProperty(prop);
+                            constexpr int nelems = 3;
+
+                            for(int e = 0; e < nelems; e++) {
+                                buffer << vector_ptr(i, e);
+                                vector_ptr(i, e) = vector_ptr(nlocal, e);
+                            }
+                        } else if(prop_type == Prop_Matrix) {
+                            auto matrix_ptr = ps->getAsMatrixProperty(prop);
+                            constexpr int nelems = 9;
+
+                            for(int e = 0; e < nelems; e++) {
+                                buffer << matrix_ptr(i, e);
+                                matrix_ptr(i, e) = matrix_ptr(nlocal, e);
+                            }
+                        } else if(prop_type == Prop_Quaternion) {
+                            auto quat_ptr = ps->getAsQuaternionProperty(prop);
+                            constexpr int nelems = 4;
+
+                            for(int e = 0; e < nelems; e++) {
+                                buffer << quat_ptr(i, e);
+                                quat_ptr(i, e) = quat_ptr(nlocal, e);
+                            }
+                        } else if(prop_type == Prop_Integer) {
+                            auto int_ptr = ps->getAsIntegerProperty(prop);
+                            buffer << int_ptr(i);
+                            int_ptr(i) = int_ptr(nlocal);
+                        } else if(prop_type == Prop_Real) {
+                            auto float_ptr = ps->getAsFloatProperty(prop);
+                            buffer << float_ptr(i);
+                            float_ptr(i) = float_ptr(nlocal);
+                        } else {
+                            std::cerr << "serializeImpl(): Invalid property type!" << std::endl;
+                            return;
                         }
-                    } else if(prop_type == Prop_Matrix) {
-                        auto matrix_ptr = ps->getAsMatrixProperty(prop);
-                        constexpr int nelems = 9;
-
-                        for(int e = 0; e < nelems; e++) {
-                            buffer << matrix_ptr(i, e);
-                            matrix_ptr(i, e) = matrix_ptr(nlocal, e);
-                        }
-                    } else if(prop_type == Prop_Quaternion) {
-                        auto quat_ptr = ps->getAsQuaternionProperty(prop);
-                        constexpr int nelems = 4;
-
-                        for(int e = 0; e < nelems; e++) {
-                            buffer << quat_ptr(i, e);
-                            quat_ptr(i, e) = quat_ptr(nlocal, e);
-                        }
-                    } else if(prop_type == Prop_Integer) {
-                        auto int_ptr = ps->getAsIntegerProperty(prop);
-                        buffer << int_ptr(i);
-                        int_ptr(i) = int_ptr(nlocal);
-                    } else if(prop_type == Prop_Real) {
-                        auto float_ptr = ps->getAsFloatProperty(prop);
-                        buffer << float_ptr(i);
-                        float_ptr(i) = float_ptr(nlocal);
-                    } else {
-                        std::cerr << "serializeImpl(): Invalid property type!" << std::endl;
-                        return 0;
                     }
                 }
-
-                // TODO: serialize contact history data as well
-                nserialized++;
             }
+
+            // TODO: serialize contact history data as well
+            nserialized++;
         }
 
         ps->setTrackedVariableAsInteger("nlocal", nlocal);
@@ -178,45 +185,46 @@ public:
         // md_resize_recv_buffer_capacity((int) nparticles);
 
         for(int i = 0; i < nrecv; ++i) {
-            for(auto &p: ps->getNonVolatileProperties()) {
-                auto prop = ps->getProperty(p_id);
-                auto prop_type = prop.getType();
+            for(auto &prop: ps->getProperties()) {
+                if(!prop.isVolatile()) {
+                    auto prop_type = prop.getType();
 
-                if(prop_type == Prop_Vector) {
-                    auto vector_ptr = ps->getAsVectorProperty(prop);
-                    constexpr int nelems = 3;
+                    if(prop_type == Prop_Vector) {
+                        auto vector_ptr = ps->getAsVectorProperty(prop);
+                        constexpr int nelems = 3;
 
-                    for(int e = 0; e < nelems; e++) {
-                        buffer >> real_tmp;
-                        vector_ptr(nlocal + i, e) = real_tmp;
-                    }
-                } else if(prop_type == Prop_Matrix) {
-                    auto matrix_ptr = ps->getAsMatrixProperty(prop);
-                    constexpr int nelems = 9;
+                        for(int e = 0; e < nelems; e++) {
+                            buffer >> real_tmp;
+                            vector_ptr(nlocal + i, e) = real_tmp;
+                        }
+                    } else if(prop_type == Prop_Matrix) {
+                        auto matrix_ptr = ps->getAsMatrixProperty(prop);
+                        constexpr int nelems = 9;
 
-                    for(int e = 0; e < nelems; e++) {
-                        buffer >> real_tmp;
-                        matrix_ptr(nlocal + i, e) = real_tmp;
-                    }
-                } else if(prop_type == Prop_Quaternion) {
-                    auto quat_ptr = ps->getAsQuaternionProperty(prop);
-                    constexpr int nelems = 4;
+                        for(int e = 0; e < nelems; e++) {
+                            buffer >> real_tmp;
+                            matrix_ptr(nlocal + i, e) = real_tmp;
+                        }
+                    } else if(prop_type == Prop_Quaternion) {
+                        auto quat_ptr = ps->getAsQuaternionProperty(prop);
+                        constexpr int nelems = 4;
 
-                    for(int e = 0; e < nelems; e++) {
+                        for(int e = 0; e < nelems; e++) {
+                            buffer >> real_tmp;
+                            quat_ptr(nlocal + i, e) = real_tmp;
+                        }
+                     } else if(prop_type == Prop_Integer) {
+                        auto int_ptr = ps->getAsIntegerProperty(prop);
+                        buffer >> int_tmp;
+                        int_ptr(nlocal + i) = int_tmp;
+                    } else if(prop_type == Prop_Real) {
+                        auto float_ptr = ps->getAsFloatProperty(prop);
                         buffer >> real_tmp;
-                        quat_ptr(nlocal + i, e) = real_tmp;
+                        float_ptr(nlocal + i) = real_tmp;
+                    } else {
+                        std::cerr << "deserializeImpl(): Invalid property type!" << std::endl;
+                        return;
                     }
-                 } else if(prop_type == Prop_Integer) {
-                    auto int_ptr = ps->getAsIntegerProperty(prop);
-                    buffer >> int_tmp;
-                    int_ptr(nlocal + i) = int_tmp;
-                } else if(prop_type == Prop_Real) {
-                    auto float_ptr = ps->getAsFloatProperty(prop);
-                    buffer >> real_tmp;
-                    float_ptr(nlocal + i) = real_tmp;
-                } else {
-                    std::cerr << "deserializeImpl(): Invalid property type!" << std::endl;
-                    return 0;
                 }
             }
         }
diff --git a/runtime/domain/block_forest.cpp b/runtime/domain/block_forest.cpp
index 9718b95ebc9711a56fdd68766119a67d964e80a8..7f0bb76fb35dbdcf5735c3fba10bfee595f3fd1b 100644
--- a/runtime/domain/block_forest.cpp
+++ b/runtime/domain/block_forest.cpp
@@ -32,7 +32,7 @@ void BlockForest::updateNeighborhood() {
 
     for(auto& iblock: *forest) {
         auto block = static_cast<walberla::blockforest::Block *>(&iblock);
-        auto& block_info = info[block->getId()];
+        auto& block_info = (*info)[block->getId()];
 
         if(block_info.computationalWeight > 0) {
             for(uint neigh = 0; neigh < block->getNeighborhoodSize(); ++neigh) {
@@ -41,7 +41,7 @@ void BlockForest::updateNeighborhood() {
                 if(neighbor_rank != me) {
                     const walberla::BlockID& neighbor_block = block->getNeighborId(neigh);
                     walberla::math::AABB neighbor_aabb = block->getNeighborAABB(neigh);
-                    auto neighbor_info = info[neighbor_block];
+                    auto neighbor_info = (*info)[neighbor_block];
                     auto begin = blocks_pushed[neighbor_rank].begin();
                     auto end = blocks_pushed[neighbor_rank].end();
 
@@ -94,11 +94,11 @@ void BlockForest::copyRuntimeArray(const std::string& name, void *dest, const in
 void BlockForest::updateWeights() {
     walberla::mpi::BufferSystem bs(walberla::mpi::MPIManager::instance()->comm(), 756);
 
-    info.clear();
+    info->clear();
     for(auto& iblock: *forest) {
         auto block = static_cast<walberla::blockforest::Block *>(&iblock);
         auto aabb = block->getAABB();
-        auto& block_info = info[block->getId()];
+        auto& block_info = (*info)[block->getId()];
 
         pairs::compute_boundary_weights(
             this->ps,
@@ -108,7 +108,7 @@ void BlockForest::updateWeights() {
         for(int branch = 0; branch < 8; ++branch) {
             const auto b_id = walberla::BlockID(block->getId(), branch);
             const auto b_aabb = forest->getAABBFromBlockId(b_id);
-            auto& b_info = info[b_id];
+            auto& b_info = (*info)[b_id];
 
             pairs::compute_boundary_weights(
                 this->ps,
@@ -119,7 +119,7 @@ void BlockForest::updateWeights() {
 
     for(auto& iblock: *forest) {
         auto block = static_cast<walberla::blockforest::Block *>(&iblock);
-        auto& block_info = info[block->getId()];
+        auto& block_info = (*info)[block->getId()];
 
         for(int neigh = 0; neigh < block->getNeighborhoodSize(); ++neigh) {
             bs.sendBuffer(block->getNeighborProcess(neigh)) <<
@@ -128,7 +128,7 @@ void BlockForest::updateWeights() {
 
         for(int branch = 0; branch < 8; ++branch) {
             const auto b_id = walberla::BlockID(block->getId(), branch);
-            auto& b_info = info[b_id];
+            auto& b_info = (*info)[b_id];
 
             for(int neigh = 0; neigh < block->getNeighborhoodSize(); ++neigh) {
                 bs.sendBuffer(block->getNeighborProcess(neigh)) <<
@@ -144,7 +144,7 @@ void BlockForest::updateWeights() {
         while(!recv.buffer().isEmpty()) {
             walberla::blockforest::InfoCollectionPair val;
             recv.buffer() >> val;
-            info.insert(val);
+            info->insert(val);
         }
     }
 }
@@ -273,30 +273,40 @@ void BlockForest::initializeWorkloadBalancer() {
     forest->setRefreshMinTargetLevelDeterminationFunction(
         walberla::blockforest::MinMaxLevelDetermination(info, regridMin, regridMax));
 
-    for_each(algorithm.begin(), algorithm.end(), [](char& c) { c = (char) ::tolower(c); });
+    std::transform(algorithm.begin(), algorithm.end(), algorithm.begin(),
+        [](unsigned char c) { return std::tolower(c); });
 
     if(algorithm == "morton") {
-        forest->setRefreshPhantomBlockDataAssignmentFunction(walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
-        forest->setRefreshPhantomBlockDataPackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
-        forest->setRefreshPhantomBlockDataUnpackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataAssignmentFunction(
+            walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
+        forest->setRefreshPhantomBlockDataPackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataUnpackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
 
         auto prepFunc = walberla::blockforest::DynamicCurveBalance<walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeight>(false, true, false);
-
         prepFunc.setMaxBlocksPerProcess(maxBlocksPerProcess);
         forest->setRefreshPhantomBlockMigrationPreparationFunction(prepFunc);
+
     } else if(algorithm == "hilbert") {
-        forest->setRefreshPhantomBlockDataAssignmentFunction(walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
-        forest->setRefreshPhantomBlockDataPackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
-        forest->setRefreshPhantomBlockDataUnpackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataAssignmentFunction(
+            walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
+        forest->setRefreshPhantomBlockDataPackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataUnpackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
 
         auto prepFunc = walberla::blockforest::DynamicCurveBalance<walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeight>(true, true, false);
-
         prepFunc.setMaxBlocksPerProcess(maxBlocksPerProcess);
         forest->setRefreshPhantomBlockMigrationPreparationFunction(prepFunc);
+
     } else if(algorithm == "metis") {
-        forest->setRefreshPhantomBlockDataAssignmentFunction(walberla::blockforest::MetisAssignmentFunctor(info, baseWeight));
-        forest->setRefreshPhantomBlockDataPackFunction(walberla::blockforest::MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
-        forest->setRefreshPhantomBlockDataUnpackFunction(walberla::blockforest::MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataAssignmentFunction(
+            walberla::blockforest::MetisAssignmentFunctor(info, baseWeight));
+        forest->setRefreshPhantomBlockDataPackFunction(
+            walberla::blockforest::MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataUnpackFunction(
+            walberla::blockforest::MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
 
         auto alg = walberla::blockforest::DynamicParMetis::stringToAlgorithm(metisAlgorithm);
         auto vWeight = walberla::blockforest::DynamicParMetis::stringToWeightsToUse(metisWeightsToUse);
@@ -305,13 +315,16 @@ void BlockForest::initializeWorkloadBalancer() {
 
         prepFunc.setipc2redist(metisipc2redist);
         forest->setRefreshPhantomBlockMigrationPreparationFunction(prepFunc);
+
     } else if(algorithm == "diffusive") {
-        forest->setRefreshPhantomBlockDataAssignmentFunction(walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
-        forest->setRefreshPhantomBlockDataPackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
-        forest->setRefreshPhantomBlockDataUnpackFunction(walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataAssignmentFunction(
+            walberla::blockforest::WeightAssignmentFunctor(info, baseWeight));
+        forest->setRefreshPhantomBlockDataPackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
+        forest->setRefreshPhantomBlockDataUnpackFunction(
+            walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor());
 
         auto prepFunc = walberla::blockforest::DynamicDiffusionBalance<walberla::blockforest::WeightAssignmentFunctor::PhantomBlockWeight>(1, 1, false);
-
         forest->setRefreshPhantomBlockMigrationPreparationFunction(prepFunc);
     }
 
diff --git a/runtime/domain/block_forest.hpp b/runtime/domain/block_forest.hpp
index 2c07fc24d7ec6cf5374669ba7b80486696b7a795..80d63d59f40d9e75389c9d650fe152e19db9cb79 100644
--- a/runtime/domain/block_forest.hpp
+++ b/runtime/domain/block_forest.hpp
@@ -26,7 +26,7 @@ class BlockForest : public DomainPartitioner {
 private:
     PairsSimulation *ps;
     std::shared_ptr<walberla::BlockForest> forest;
-    walberla::blockforest::InfoCollection info;
+    std::shared_ptr<walberla::blockforest::InfoCollection> info;
     std::map<int, std::vector<walberla::math::AABB>> neighborhood;
     std::map<int, std::vector<walberla::BlockID>> blocks_pushed;
     std::vector<int> ranks;
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 8edad7f772b4433aabee47bd931c47ee00a0291d..9af0f9737707b239e11f21363ac656511fab2c2f 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -28,7 +28,6 @@ namespace pairs {
 
 class PairsSimulation {
 private:
-    //Regular6DStencil *dom_part;
     DomainPartitioner *dom_part;
     DomainPartitioners dom_part_type;
     std::vector<Property> properties;
@@ -147,6 +146,7 @@ public:
     void copyArraySliceToHost(Array &array, action_t action, size_t offset, size_t size);
 
     // Properties
+    std::vector<Property> &getProperties() { return properties; };
     Property &getProperty(property_t id);
     Property &getPropertyByName(std::string name);
     void addProperty(Property prop);
@@ -154,11 +154,11 @@ public:
     template<typename T_ptr>
     void addProperty(
         property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t,
-        PropertyType type, layout_t layout, size_t sx, size_t sy = 1);
+        PropertyType type, layout_t layout, int vol, size_t sx, size_t sy = 1);
 
     template<typename T_ptr> void addProperty(
         property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr,
-        PropertyType type, layout_t layout, size_t sx, size_t sy = 1);
+        PropertyType type, layout_t layout, int vol, size_t sx, size_t sy = 1);
 
     template<typename T_ptr>
     void reallocProperty(property_t id, T_ptr **h_ptr, std::nullptr_t, size_t sx = 1, size_t sy = 1);
@@ -405,19 +405,21 @@ void PairsSimulation::reallocArray(array_t id, T_ptr **h_ptr, T_ptr **d_ptr, siz
 
 template<typename T_ptr>
 void PairsSimulation::addProperty(
-    property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, PropertyType type, layout_t layout, size_t sx, size_t sy) {
+    property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t,
+    PropertyType type, layout_t layout, int vol, size_t sx, size_t sy) {
 
     size_t size = sx * sy * sizeof(T_ptr);
     PAIRS_ASSERT(size > 0);
 
     *h_ptr = (T_ptr *) pairs::host_alloc(size);
     PAIRS_ASSERT(*h_ptr != nullptr);
-    addProperty(Property(id, name, *h_ptr, nullptr, type, layout, sx, sy));
+    addProperty(Property(id, name, *h_ptr, nullptr, type, layout, vol, sx, sy));
 }
 
 template<typename T_ptr>
 void PairsSimulation::addProperty(
-    property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, PropertyType type, layout_t layout, size_t sx, size_t sy) {
+    property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr,
+    PropertyType type, layout_t layout, int vol, size_t sx, size_t sy) {
 
     size_t size = sx * sy * sizeof(T_ptr);
     PAIRS_ASSERT(size > 0);
@@ -425,7 +427,7 @@ void PairsSimulation::addProperty(
     *h_ptr = (T_ptr *) pairs::host_alloc(size);
     *d_ptr = (T_ptr *) pairs::device_alloc(size);
     PAIRS_ASSERT(*h_ptr != nullptr && *d_ptr != nullptr);
-    addProperty(Property(id, name, *h_ptr, *d_ptr, type, layout, sx, sy));
+    addProperty(Property(id, name, *h_ptr, *d_ptr, type, layout, vol, sx, sy));
 }
 
 template<typename T_ptr>
diff --git a/runtime/property.hpp b/runtime/property.hpp
index 741594d1745edf923f63bf8eae422aa11218f41a..84a2bc38151ee8ada0a0ebcf2283f8d5d5303bc5 100644
--- a/runtime/property.hpp
+++ b/runtime/property.hpp
@@ -12,15 +12,20 @@ protected:
     PropertyType type;
     layout_t layout;
     size_t sx, sy;
+    int vol;
 
 public:
-    Property(property_t id_, std::string name_, void *h_ptr_, void *d_ptr_, PropertyType type_, layout_t layout_, size_t sx_, size_t sy_=1) :
+    Property(
+        property_t id_, std::string name_, void *h_ptr_, void *d_ptr_,
+        PropertyType type_, layout_t layout_, int vol_, size_t sx_, size_t sy_=1) :
+
         id(id_),
         name(name_),
         h_ptr(h_ptr_),
         d_ptr(d_ptr_),
         type(type_),
         layout(layout_),
+        vol(vol_),
         sx(sx_), sy(sy_) {
 
         PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
@@ -35,6 +40,7 @@ public:
     size_t getTotalSize() { return sx * sy * getPrimitiveTypeSize(); };
     PropertyType getType() { return type; }
     layout_t getLayout() { return layout; }
+    int isVolatile() { return vol != 0; }
     size_t getPrimitiveTypeSize() {
         return  (type == Prop_Integer) ? sizeof(int) :
                 (type == Prop_Real) ? sizeof(real_t) :
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 6e1b310e7a2c83de0c266ac58115900208ee38e7..226af7f6c99d81837c2dc5656be23766f8623cfe 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -698,6 +698,7 @@ class CGen:
             assert ptype != "Prop_Invalid", "Invalid property type!"
 
             playout = Layouts.c_keyword(p.layout())
+            vol = 1 if p.is_volatile() else 0
             sizes = ", ".join([str(self.generate_expression(ScalarOp.inline(size))) for size in ast_node.sizes()])
 
             if self.target.is_gpu() and p.device_flag:
@@ -706,7 +707,7 @@ class CGen:
             else:
                 self.print(f"{tkw} *{ptr};")
 
-            self.print(f"pairs->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
+            self.print(f"pairs->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {vol}, {sizes});")
 
         if isinstance(ast_node, RegisterContactProperty):
             p = ast_node.property()
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index f10f9424b82f721b3564f6fdc9899f94aca4bc7e..b13d6f48e8face187f328130a95a14eb1ff0717d 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -83,6 +83,9 @@ class Property(ASTNode):
     def default(self):
         return self.default_value
 
+    def is_volatile(self):
+        return self.volatile
+
     def ndims(self):
         return 1 if Types.is_scalar(self.prop_type) else 2