diff --git a/src/mesa_pd/kernel/cnt/IntegratedVDWContact.h b/src/mesa_pd/kernel/cnt/IntegratedVDWContact.h
new file mode 100644
index 0000000000000000000000000000000000000000..844887f15c5b472b12830f0e933709f4e69681b1
--- /dev/null
+++ b/src/mesa_pd/kernel/cnt/IntegratedVDWContact.h
@@ -0,0 +1,145 @@
+//======================================================================================================================
+//
+// 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
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+#include <core/math/Constants.h>
+#include <core/logging/Logging.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+namespace cnt {
+
+/**
+ * vdW Contact with integration
+ */
+class IntegratedVDWContact
+{
+public:
+ template<typename Accessor>
+ void operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor &ac);
+
+ static constexpr real_t R_CNT = 6.78_r; ///< CNT radius
+ static constexpr real_t T = 2_r * R_CNT; ///< Height of a cylindrical segment
+ static constexpr real_t eps = 0.07124_r;
+ static constexpr real_t A = 0.0223_r;
+ static constexpr real_t B = 1.31_r;
+ static constexpr real_t alpha = 9.5_r;
+ static constexpr real_t beta = 4.0_r;
+
+ // vdW adhesion + linear repulsion potential.
+ const real_t r0 = R_CNT * (2_r + std::pow((alpha * A / (beta * B)), 1_r / (alpha - beta)));
+ const real_t u0 = 4_r * eps * (A / std::pow(r0 / R_CNT - 2_r, (alpha)) - B / std::pow(r0 / R_CNT - 2_r, (beta)));
+
+ auto getLastEnergy() const {return energy_;}
+private:
+ real_t energy_; ///< total potential
+};
+
+
+template<typename Accessor>
+inline
+void IntegratedVDWContact::operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor &ac)
+{
+ constexpr real_t K_n = 200.0_r; // Good for fabrics modeling
+
+ // particle centers
+ Vec3 O1 = ac.getPosition(p_idx1);
+ Vec3 O2 = ac.getPosition(p_idx2);
+
+ // axial vectors
+ Vec3 b1 = ac.getRotation(p_idx1).getMatrix() * Vec3(1_r, 0_r, 0_r);
+ Vec3 b2 = ac.getRotation(p_idx2).getMatrix() * Vec3(1_r, 0_r, 0_r);
+
+ energy_ = 0_r;
+ Vec3 force12(0); // Force 12
+ Vec3 force21(0); // Force 21
+ Vec3 moment12(0); // Total torque 12
+ Vec3 moment21(0); // Total torque 21
+
+ constexpr int Np = 5; // Number of integration points over each axis
+ constexpr real_t Npinv = 1.0_r / real_t(Np);
+ for (int i = 0; i < Np; ++i) // integral dl1
+ {
+ for (int j = 0; j < Np; ++j) // integral dl2
+ {
+ // Levers
+ Vec3 l1 = (-0.5_r * T + (0.5_r + real_t(i)) * T * Npinv) * b1;
+ Vec3 l2 = (-0.5_r * T + (0.5_r + real_t(j)) * T * Npinv) * b2;
+
+ /// radius vector between dl1 and dl2
+ Vec3 R12 = (O2 + l2) - (O1 + l1);
+
+ real_t r12 = R12.length();
+ Vec3 n12 = R12 * (1_r / r12);
+
+ real_t dU = 0_r;
+ Vec3 dforce12(0);
+
+ if (r12 < r0)
+ {
+ // elastic interaction
+ dU = u0 + K_n * (r12 - r0) * (r12 - r0) * Npinv * Npinv;
+ dforce12 = n12 * K_n * (r12 - r0) * Npinv * Npinv;
+ } else
+ {
+ // vdW interaction
+ const real_t normDistance = r12 / R_CNT - 2_r;
+ const real_t powAlpha = std::pow(normDistance, alpha);
+ const real_t powBeta = std::pow(normDistance, beta);
+ dU = 4_r * eps * (A / powAlpha - B / powBeta) * Npinv * Npinv;
+ dforce12 = n12 * 4_r * eps / R_CNT * Npinv * Npinv *
+ (-(alpha * A) / (powAlpha * normDistance) + (beta * B) / (powBeta * normDistance));
+ }
+
+ Vec3 dmoment12 = l2 % dforce12;
+ Vec3 dmoment21 = l1 % (-dforce12);
+
+ energy_ += dU;
+ force12 += dforce12;
+ force21 -= dforce12;
+ moment12 += dmoment12;
+ moment21 += dmoment21;
+ }
+ }
+
+ addForceAtomic(p_idx1, ac, force12);
+ addForceAtomic(p_idx2, ac, force21);
+ addTorqueAtomic(p_idx1, ac, -moment21);
+ addTorqueAtomic(p_idx2, ac, -moment12);
+}
+
+} //namespace cnt
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/src/mesa_pd/kernel/cnt/Parameters.h b/src/mesa_pd/kernel/cnt/Parameters.h
index 63f256a04dca173007766af03f3b70256047824f..c265088bab02da9837599fbf333032b795bf2b0d 100644
--- a/src/mesa_pd/kernel/cnt/Parameters.h
+++ b/src/mesa_pd/kernel/cnt/Parameters.h
@@ -79,11 +79,11 @@ constexpr auto M_C = 12.011_r;
/// Mass of the repetative cell in AMU
constexpr auto mass_T = ro * T * M_C * 104.397_r;
-// Volume of a capsule
-//double vol = (4./3.) * PII * (R_CNT * R_CNT * R_CNT) + PII * R_CNT * R_CNT * T;
+/// Volume of a capsule
+double vol_capsule = (4_r/3_r) * math::pi * (R_CNT * R_CNT * R_CNT) + math::pi * R_CNT * R_CNT * T;
-// Density of a capsule
-//double dens = mass_T / vol;
+/// Density of a capsule
+double dens_capsule = mass_T / vol_capsule;
/// V-bond parameter
constexpr auto knorm = (En * 0.006242_r) / T;
@@ -97,12 +97,12 @@ constexpr auto outer_radius = CutoffFactor * R_CNT;
constexpr auto inner_radius = 1.5811_r * R_CNT; // Sphere
/// volume of a sphere
-constexpr auto vol = (4_r/3_r) * math::pi * (inner_radius * inner_radius * inner_radius);
+constexpr auto vol_sphere = (4_r/3_r) * math::pi * (inner_radius * inner_radius * inner_radius);
/// density of a sphere
-constexpr auto dens = mass_T / vol;
+constexpr auto dens_sphere = mass_T / vol_sphere;
} //namespace cnt
} //namespace kernel
} //namespace mesa_pd
-} //namespace walberla
\ No newline at end of file
+} //namespace walberla
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 983cbca6e994be63a8cc4a275d6e8457b501a24a..4e7985268047a837470ea5ed0b2a4793b1785873 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -84,6 +84,9 @@ waLBerla_execute_test( NAME MESA_PD_DropTestGeneral )
waLBerla_compile_test( NAME MESA_PD_Kernel_ClearNextNeighborSync FILES kernel/ClearNextNeighborSync.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_ClearNextNeighborSync PROCESSES 2 )
+waLBerla_compile_test( NAME MESA_PD_Kernel_CNT_IntegratedVDWContact FILES kernel/cnt/IntegratedVDWContact.test.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_CNT_IntegratedVDWContact )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_CNT_IsotropicVDWContact FILES kernel/cnt/IsotropicVDWContact.test.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_CNT_IsotropicVDWContact )
diff --git a/tests/mesa_pd/kernel/cnt/IntegratedVDWContact.test.cpp b/tests/mesa_pd/kernel/cnt/IntegratedVDWContact.test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..05bb96c988d9785bce0b12de23ffb327bb762893
--- /dev/null
+++ b/tests/mesa_pd/kernel/cnt/IntegratedVDWContact.test.cpp
@@ -0,0 +1,135 @@
+//======================================================================================================================
+//
+// 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
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "SphericalSegmentAccessor.h"
+
+#include "mesa_pd/data/Flags.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/cnt/IntegratedVDWContact.h"
+#include "mesa_pd/kernel/cnt/ViscousDamping.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "core/Environment.h"
+#include "core/math/Constants.h"
+#include "vtk/VTKOutput.h"
+
+namespace walberla {
+using namespace walberla::mesa_pd;
+
+int main(int argc, char **argv)
+{
+ Environment env(argc, argv);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ if (std::is_same<walberla::real_t, float>::value)
+ {
+ WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+ return EXIT_SUCCESS;
+ }
+
+ logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+ WALBERLA_LOG_INFO_ON_ROOT("loading configuration parameters");
+ constexpr auto numSimulationSteps = 20000ll;
+ constexpr auto outputInterval = 100ll;
+
+ WALBERLA_LOG_INFO_ON_ROOT("creating initial particle setup");
+ auto ps = std::make_shared<data::ParticleStorage>(10);
+ auto ac = SphericalSegmentAccessor(ps);
+
+ using namespace kernel::cnt;
+ data::Particle &&sp1 = *ps->create();
+ sp1.setPosition(Vec3(0_r, 0_r, 0_r));
+ sp1.setSegmentID(1);
+ sp1.setClusterID(1);
+
+ data::Particle &&sp2 = *ps->create();
+ sp2.setPosition(Vec3(20_r, 20_r, 20_r));
+ sp2.setSegmentID(2);
+ sp2.setClusterID(2);
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up VTK output");
+ auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+ vtkOutput->addOutput<data::SelectParticlePosition>("position");
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput,
+ "cnt",
+ 1,
+ "vtk_integrated",
+ "particles",
+ false,
+ false);
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up interaction models");
+ kernel::cnt::IntegratedVDWContact vdW_integrated;
+ kernel::cnt::ViscousDamping viscous_damping(0.1_r * 1052.0_r, 0.1_r * 1052.0_r);
+ kernel::VelocityVerletPreForceUpdate vv_pre(kernel::cnt::dT);
+ kernel::VelocityVerletPostForceUpdate vv_post(kernel::cnt::dT);
+
+ WALBERLA_LOG_INFO_ON_ROOT("running simulation");
+
+ real_t U = 0_r;
+ for (auto i = 0; i < numSimulationSteps; ++i)
+ {
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_pre,
+ ac);
+
+ U = 0_r;
+ ps->forEachParticlePairHalf(false,
+ kernel::SelectAll(),
+ ac,
+ [&](size_t p_idx1, size_t p_idx2)
+ {
+ vdW_integrated(p_idx1, p_idx2, ac);
+ U += vdW_integrated.getLastEnergy();
+ viscous_damping(p_idx1, p_idx2, ac);
+ });
+
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_post,
+ ac);
+
+ if( i % outputInterval == 0 )
+ {
+// vtkWriter->write();
+// WALBERLA_LOG_DEVEL(i << " : " << U);
+ }
+ }
+
+ WALBERLA_CHECK_LESS(U, -8_r)
+
+ return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char *argv[])
+{
+ return walberla::main(argc, argv);
+}