diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index bb3cb83c6279ccbd2e356747199045996afcebf6..bd8baa406fc93a5d2ff4e7869fa79e5d278c2181 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -15,6 +15,7 @@ if __name__ == '__main__':
     args = parser.parse_args()
 
     mpd = Module(args.path)
+    mpd.enable_openmp(False)
     ps = mpd.add(data.ParticleStorage())
     ps.set_shapes('Sphere', 'HalfSpace', 'CylindricalBoundary', 'Box', 'Ellipsoid')
     ps.add_property("position", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
diff --git a/python/mesa_pd/Module.py b/python/mesa_pd/Module.py
index 64c307c4674d04524377e1f05df4fb5024626a11..57d989ffa2df7e57cd794fbaf8b11a9593b7a080 100644
--- a/python/mesa_pd/Module.py
+++ b/python/mesa_pd/Module.py
@@ -17,11 +17,14 @@ class Module:
            name of the generated module
         """
 
-        self.context = {}
-        self.context['path'] = Path(path).resolve()
-        self.context['name'] = module_name
-        self.context['module_path'] = self.context['path'] / 'src' / self.context['name']
-        self.context['test_path']   = self.context['path'] / 'tests' / self.context['name']
+        path = Path(path).resolve()
+        self.context = {
+            'path': path,
+            'name': module_name,
+            'module_path': path / 'src' / module_name,
+            'test_path': path / 'tests' / module_name,
+            'enableOpenMP': False
+        }
 
         self.components = []
 
@@ -29,6 +32,9 @@ class Module:
         self.components.append(component)
         return component
 
+    def enable_openmp(self, enabled):
+        self.context['enableOpenMP'] = enabled
+
     def rename(self):
         for filename in (f for f in self.context['module_path'].glob('**/*') if f.is_file()):
             filedata = None
diff --git a/python/mesa_pd/data/ParticleStorage.py b/python/mesa_pd/data/ParticleStorage.py
index 65193bab26d9f8f8437d2c54ca5b09badc8dea92..76e840dff6eaca7da936110e4cd7668a626a6802 100644
--- a/python/mesa_pd/data/ParticleStorage.py
+++ b/python/mesa_pd/data/ParticleStorage.py
@@ -100,5 +100,6 @@ class ParticleStorage():
     def generate(self, module):
         ctx = {'module': module, **self.context}
 
-        generate_file(module['module_path'], 'data/ParticleStorage.templ.h', ctx, filename='data/ParticleStorage.h')
-        generate_file(module['module_path'], 'data/ParticleAccessor.templ.h', ctx, filename='data/ParticleAccessor.h')
+        generate_file(module['module_path'], 'data/ParticleStorage.templ.h', ctx)
+        generate_file(module['module_path'], 'data/ParticleAccessor.templ.h', ctx)
+        generate_file(module['module_path'], 'common/ParticleFunctions.templ.h', ctx)
diff --git a/python/mesa_pd/templates/common/ParticleFunctions.templ.h b/python/mesa_pd/templates/common/ParticleFunctions.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ed2d53db7188722ef1ad457105c9e32202576bd
--- /dev/null
+++ b/python/mesa_pd/templates/common/ParticleFunctions.templ.h
@@ -0,0 +1,148 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ParticleFunctions.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+//  THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/domain/IDomain.h>
+
+#include <core/logging/Logging.h>
+#include <core/mpi/MPIManager.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+/**
+ * Returns the "surface" velocity at a certain point given in world frame coordinates.
+ */
+template <typename Accessor>
+inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& wf_pt)
+{
+   return ac.getLinearVelocity(p_idx) + cross(ac.getAngularVelocity(p_idx), ( wf_pt - ac.getPosition(p_idx) ));
+}
+
+/**
+ * Force is applied at the center of mass.
+ */
+template <typename Accessor>
+inline void addForceAtomic(const size_t p_idx, Accessor& ac, const Vec3& f)
+{
+   // Increasing the force and torque on this particle
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[0]  += f[0];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[1]  += f[1];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[2]  += f[2];
+}
+
+template <typename Accessor>
+inline void addForceAtWFPosAtomic(const size_t p_idx, Accessor& ac, const Vec3& f, const Vec3& wf_pt)
+{
+   // Increasing the force and torque on this particle
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[0]  += f[0];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[1]  += f[1];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getForceRef(p_idx)[2]  += f[2];
+
+   const auto t = cross(( wf_pt - ac.getPosition(p_idx) ), f);
+
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[0] += t[0];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[1] += t[1];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[2] += t[2];
+}
+
+/**
+ * Torque is directly applied on the particle.
+ */
+template <typename Accessor>
+inline void addTorqueAtomic(const size_t p_idx, Accessor& ac, const Vec3& t)
+{
+   // Increasing the torque on this particle
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[0]  += t[0];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[1]  += t[1];
+   {%- if module.enableOpenMP %}
+   #ifdef _OPENMP
+   #pragma omp atomic
+   #endif
+   {%- endif %};
+   ac.getTorqueRef(p_idx)[2]  += t[2];
+}
+
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd/templates/data/ContactStorage.templ.h b/python/mesa_pd/templates/data/ContactStorage.templ.h
index f0af2d1fa8c28264c0b3ee719631d4a6a2cb67eb..3731188bb421575cd41242f71186c891b5bd365f 100644
--- a/python/mesa_pd/templates/data/ContactStorage.templ.h
+++ b/python/mesa_pd/templates/data/ContactStorage.templ.h
@@ -369,9 +369,9 @@ inline void ContactStorage::forEachContact(const bool openmp, const Selector& se
 {
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for schedule(static) firstprivate(selector,func) if (openmp)
+   {%- endif %}
    for (int64_t i = 0; i < int64_c(len); ++i)
       if (selector(uint64_c(i), acForPS)){
          func( uint64_c(i), std::forward<Args>(args)... );
diff --git a/python/mesa_pd/templates/data/LinkedCells.templ.h b/python/mesa_pd/templates/data/LinkedCells.templ.h
index dc13cfdf4efbc0f790c20360f03d9360b8fe33ab..713e2d073810c76029212d19010e7da87c344d00 100644
--- a/python/mesa_pd/templates/data/LinkedCells.templ.h
+++ b/python/mesa_pd/templates/data/LinkedCells.templ.h
@@ -171,9 +171,9 @@ void LinkedCells::clear()
 {
    const uint64_t cellsSize = cells_.size();
    //clear existing linked cells
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static)
-#endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for schedule(static)
+   {%- endif %}
    for (int64_t i = 0; i < int64_c(cellsSize); ++i)
       cells_[uint64_c(i)] = -1;
    infiniteParticles_ = -1;
@@ -185,9 +185,9 @@ inline void LinkedCells::forEachParticlePair{%- if half %}Half{%- endif %}(const
 {
    static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
    WALBERLA_UNUSED(openmp);
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static) if (openmp)
-#endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for collapse(2) schedule(static) firstprivate(selector,func) if (openmp)
+   {%- endif %}
    for (int z = 0; z < numCellsPerDim_[2]; ++z)
    {
       for (int y = 0; y < numCellsPerDim_[1]; ++y)
diff --git a/python/mesa_pd/templates/data/ParticleStorage.templ.h b/python/mesa_pd/templates/data/ParticleStorage.templ.h
index 72b309e241404bf9b7bc5e14c195ddf527122df2..7f8d00bb3b5884e17dbfc57a4cf20596dbb373fe 100644
--- a/python/mesa_pd/templates/data/ParticleStorage.templ.h
+++ b/python/mesa_pd/templates/data/ParticleStorage.templ.h
@@ -490,9 +490,9 @@ inline void ParticleStorage::forEachParticle(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for schedule(static) firstprivate(selector,func) if (openmp)
+   {%- endif %}
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       if (selector(uint64_c(i), acForPS))
@@ -512,9 +512,9 @@ inline void ParticleStorage::forEachParticlePair(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for schedule(static) firstprivate(selector,func) if (openmp)
+   {%- endif %}
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = 0; j < int64_c(len); ++j)
@@ -542,9 +542,9 @@ inline void ParticleStorage::forEachParticlePairHalf(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
+   {%- if module.enableOpenMP %}
+   #pragma omp parallel for schedule(static) firstprivate(selector,func) if (openmp)
+   {%- endif %}
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = i+1; j < int64_c(len); ++j)
diff --git a/python/mesa_pd/templates/kernel/ForceLJ.templ.h b/python/mesa_pd/templates/kernel/ForceLJ.templ.h
index c0b6234a380e9abb76b20745ee6c08d715e8bb2d..4242a0658d0fa2e8e4007f78349616deca5d0ab8 100644
--- a/python/mesa_pd/templates/kernel/ForceLJ.templ.h
+++ b/python/mesa_pd/templates/kernel/ForceLJ.templ.h
@@ -26,6 +26,7 @@
 
 #pragma once
 
+#include <mesa_pd/common/ParticleFunctions.h>
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
 
@@ -129,31 +130,9 @@ inline void ForceLJ::operator()(const size_t p_idx, const size_t np_idx, Accesso
       const real_t force = real_t(48) * sr6 * ( sr6 - real_t(0.5) ) * sr2 * getEpsilon(ac.getType(p_idx), ac.getType(np_idx));
       const Vec3 f = force * dir;
 
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[0]  += f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[1]  += f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[2]  += f[2];
-
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[0]  -= f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[1]  -= f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[2]  -= f[2];
+      // Add normal force at contact point
+      addForceAtomic( p_idx, ac, f );
+      addForceAtomic( np_idx, ac, -f );
    }
 }
 
diff --git a/src/mesa_pd/common/ParticleFunctions.h b/src/mesa_pd/common/ParticleFunctions.h
index d19045ab05c4cce40858672b5572188a337c8a5c..462e3f2e8ec428b10d56171e54f1baf97597c11c 100644
--- a/src/mesa_pd/common/ParticleFunctions.h
+++ b/src/mesa_pd/common/ParticleFunctions.h
@@ -44,51 +44,23 @@ inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& w
 template <typename Accessor>
 inline void addForceAtomic(const size_t p_idx, Accessor& ac, const Vec3& f)
 {
-   // Increasing the force and torque on this particle
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[0]  += f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[1]  += f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
+   // Increasing the force and torque on this particle;
+   ac.getForceRef(p_idx)[0]  += f[0];;
+   ac.getForceRef(p_idx)[1]  += f[1];;
    ac.getForceRef(p_idx)[2]  += f[2];
 }
 
 template <typename Accessor>
 inline void addForceAtWFPosAtomic(const size_t p_idx, Accessor& ac, const Vec3& f, const Vec3& wf_pt)
 {
-   // Increasing the force and torque on this particle
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[0]  += f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[1]  += f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
+   // Increasing the force and torque on this particle;
+   ac.getForceRef(p_idx)[0]  += f[0];;
+   ac.getForceRef(p_idx)[1]  += f[1];;
    ac.getForceRef(p_idx)[2]  += f[2];
 
-   const auto t = cross(( wf_pt - ac.getPosition(p_idx) ), f);
-
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getTorqueRef(p_idx)[0] += t[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getTorqueRef(p_idx)[1] += t[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
+   const auto t = cross(( wf_pt - ac.getPosition(p_idx) ), f);;
+   ac.getTorqueRef(p_idx)[0] += t[0];;
+   ac.getTorqueRef(p_idx)[1] += t[1];;
    ac.getTorqueRef(p_idx)[2] += t[2];
 }
 
@@ -98,21 +70,12 @@ inline void addForceAtWFPosAtomic(const size_t p_idx, Accessor& ac, const Vec3&
 template <typename Accessor>
 inline void addTorqueAtomic(const size_t p_idx, Accessor& ac, const Vec3& t)
 {
-   // Increasing the torque on this particle
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getTorqueRef(p_idx)[0]  += t[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getTorqueRef(p_idx)[1]  += t[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
+   // Increasing the torque on this particle;
+   ac.getTorqueRef(p_idx)[0]  += t[0];;
+   ac.getTorqueRef(p_idx)[1]  += t[1];;
    ac.getTorqueRef(p_idx)[2]  += t[2];
 }
 
 
 } //namespace mesa_pd
-} //namespace walberla
+} //namespace walberla
\ No newline at end of file
diff --git a/src/mesa_pd/data/ContactStorage.h b/src/mesa_pd/data/ContactStorage.h
index b257a59a0269118d5c7b65f5d737847cc40de72c..0de025bac7d673b400a0fc64d74a5d5cb724d629 100644
--- a/src/mesa_pd/data/ContactStorage.h
+++ b/src/mesa_pd/data/ContactStorage.h
@@ -601,9 +601,6 @@ inline void ContactStorage::forEachContact(const bool openmp, const Selector& se
 {
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
       if (selector(uint64_c(i), acForPS)){
          func( uint64_c(i), std::forward<Args>(args)... );
@@ -615,9 +612,6 @@ inline void ContactStorage::forEachContact(const bool openmp, const Selector& se
 {
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
       if (selector(uint64_c(i), acForPS)){
          func( uint64_c(i), std::forward<Args>(args)... );
diff --git a/src/mesa_pd/data/LinkedCells.h b/src/mesa_pd/data/LinkedCells.h
index f3d4ea54d32a05b6a87c2940d8145a850f71f358..3915558a90a696099dff12575bd055ca2def9a06 100644
--- a/src/mesa_pd/data/LinkedCells.h
+++ b/src/mesa_pd/data/LinkedCells.h
@@ -175,9 +175,6 @@ void LinkedCells::clear()
 {
    const uint64_t cellsSize = cells_.size();
    //clear existing linked cells
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static)
-#endif
    for (int64_t i = 0; i < int64_c(cellsSize); ++i)
       cells_[uint64_c(i)] = -1;
    infiniteParticles_ = -1;
@@ -187,9 +184,6 @@ inline void LinkedCells::forEachParticlePair(const bool openmp, const Selector&
 {
    static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
    WALBERLA_UNUSED(openmp);
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static) if (openmp)
-#endif
    for (int z = 0; z < numCellsPerDim_[2]; ++z)
    {
       for (int y = 0; y < numCellsPerDim_[1]; ++y)
@@ -282,9 +276,6 @@ inline void LinkedCells::forEachParticlePairHalf(const bool openmp, const Select
 {
    static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
    WALBERLA_UNUSED(openmp);
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static) if (openmp)
-#endif
    for (int z = 0; z < numCellsPerDim_[2]; ++z)
    {
       for (int y = 0; y < numCellsPerDim_[1]; ++y)
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index 8e248aca12c376fbf214481a7f9c552ffd5ccfcc..f4f8de7e5afb83145bf35d244a5b139c1580a98a 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -1045,9 +1045,6 @@ inline void ParticleStorage::forEachParticle(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       if (selector(uint64_c(i), acForPS))
@@ -1064,9 +1061,6 @@ inline void ParticleStorage::forEachParticle(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       if (selector(uint64_c(i), acForPS))
@@ -1083,9 +1077,6 @@ inline void ParticleStorage::forEachParticlePair(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = 0; j < int64_c(len); ++j)
@@ -1110,9 +1101,6 @@ inline void ParticleStorage::forEachParticlePair(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = 0; j < int64_c(len); ++j)
@@ -1137,9 +1125,6 @@ inline void ParticleStorage::forEachParticlePairHalf(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = i+1; j < int64_c(len); ++j)
@@ -1161,9 +1146,6 @@ inline void ParticleStorage::forEachParticlePairHalf(const bool openmp,
 
    WALBERLA_UNUSED(openmp);
    const uint64_t len = size();
-   #ifdef _OPENMP
-   #pragma omp parallel for schedule(static) if (openmp)
-   #endif
    for (int64_t i = 0; i < int64_c(len); ++i)
    {
       for (int64_t j = i+1; j < int64_c(len); ++j)
diff --git a/src/mesa_pd/kernel/ForceLJ.h b/src/mesa_pd/kernel/ForceLJ.h
index dd40c699a766dbc24079d272365ed4a8fced8cd3..b937372b96ffeaf09bca1f69480426050cd79cd3 100644
--- a/src/mesa_pd/kernel/ForceLJ.h
+++ b/src/mesa_pd/kernel/ForceLJ.h
@@ -26,6 +26,7 @@
 
 #pragma once
 
+#include <mesa_pd/common/ParticleFunctions.h>
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
 
@@ -138,31 +139,9 @@ inline void ForceLJ::operator()(const size_t p_idx, const size_t np_idx, Accesso
       const real_t force = real_t(48) * sr6 * ( sr6 - real_t(0.5) ) * sr2 * getEpsilon(ac.getType(p_idx), ac.getType(np_idx));
       const Vec3 f = force * dir;
 
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[0]  += f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[1]  += f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(p_idx)[2]  += f[2];
-
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[0]  -= f[0];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[1]  -= f[1];
-#ifdef _OPENMP
-#pragma omp atomic
-#endif
-   ac.getForceRef(np_idx)[2]  -= f[2];
+      // Add normal force at contact point
+      addForceAtomic( p_idx, ac, f );
+      addForceAtomic( np_idx, ac, -f );
    }
 }
 
diff --git a/tests/mesa_pd/kernel/ForceLJ.cpp b/tests/mesa_pd/kernel/ForceLJ.cpp
index 6e15f77740334e783eab1ecc12860c528514c8df..12c49b5c7aa38b3c4f85b51afa723b185b993fe5 100644
--- a/tests/mesa_pd/kernel/ForceLJ.cpp
+++ b/tests/mesa_pd/kernel/ForceLJ.cpp
@@ -102,21 +102,6 @@ int main( int argc, char ** argv )
    lj(0, 1, accessor);
    WALBERLA_CHECK_FLOAT_EQUAL( p1.getForce(), -p2.getForce() );
 
-   // thread safety test
-   auto singleForce = p1.getForce();
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static)
-#endif
-   for (int i = 0; i < 100; ++i)
-      lj(0, 1, accessor);
-
-   WALBERLA_LOG_DEVEL(p1);
-   WALBERLA_LOG_DEVEL(p2);
-   WALBERLA_CHECK_FLOAT_EQUAL( p1.getForce(), -p2.getForce() );
-   WALBERLA_CHECK_FLOAT_EQUAL_EPSILON( p1.getForce(),
-                                       real_t(101) * singleForce,
-                                       real_t(1e-13) );
-
    return EXIT_SUCCESS;
 }
 
diff --git a/tests/mesa_pd/kernel/SpringDashpot.cpp b/tests/mesa_pd/kernel/SpringDashpot.cpp
index 8bd82ee16481e95aa763d3ed467ce41df36ae03e..b9feabe924b849106638efc62994f5a93d0758eb 100644
--- a/tests/mesa_pd/kernel/SpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/SpringDashpot.cpp
@@ -107,23 +107,6 @@ int main( int argc, char ** argv )
    WALBERLA_CHECK_FLOAT_EQUAL( ps->getForce(0), -ps->getForce(1) );
    WALBERLA_CHECK_FLOAT_EQUAL( ps->getForce(0), Vec3(1,1,1).getNormalized() * ((std::sqrt(real_t(12)) - 4) * sd.getStiffness(0, 0)) );
 
-   // thread safety test
-#ifdef _OPENMP
-#pragma omp parallel for schedule(static)
-#endif
-   for (int i = 0; i < 100; ++i)
-      sd(contact.getIdx1(),
-         contact.getIdx2(),
-         ac,
-         contact.getContactPoint(),
-         contact.getContactNormal(),
-         contact.getPenetrationDepth());
-
-   WALBERLA_CHECK_FLOAT_EQUAL( ps->getForce(0), -ps->getForce(1) );
-   WALBERLA_CHECK_FLOAT_EQUAL_EPSILON( ps->getForce(0),
-                                       real_t(101) * Vec3(1,1,1).getNormalized() * ((std::sqrt(real_t(12)) - 4) * sd.getStiffness(0, 0)),
-                                       real_t(1e-6) );
-
    auto cor  = real_t(0.87);
    auto ct   = real_t(0.17);
    auto meff = real_t(0.65);