Code maintenance (espressomd#4562)

Description of changes:
- skip statistical tests in the macOS CI job to reduce runtime from 1h30 to 1h
- cleanup header includes and feature guards (fixes espressomd#4560)

.github/actions/build_and_check/action.yml

@@ -9,10 +9,6 @@ inputs:
     description: 'Whether to build with undefined behavior sanitizer'
     required: true
     default: 'false'
-  check_skip_long:  # id of input
-    description: 'Whether to skip long python tests'
-    required: true
-    default: 'false'
 runs:
   using: "composite"
   steps:
@@ -22,7 +18,7 @@ runs:
        pip3 install -c requirements.txt numpy cython h5py scipy
       shell: bash
     - run: |
-        export myconfig=maxset with_cuda=false test_timeout=800 with_asan=${{ inputs.asan }} with_ubsan=${{ inputs.ubsan }} check_skip_long=${{ inputs.check_skip_long }}
+        export myconfig=maxset with_cuda=false test_timeout=800 with_asan=${{ inputs.asan }} with_ubsan=${{ inputs.ubsan }} check_skip_long=true
         bash maintainer/CI/build_cmake.sh
       shell: bash
       # This is a workaround for the unfortunate interaction of MacOS and OpenMPI 4

.github/workflows/push_pull.yml

@@ -25,7 +25,6 @@ jobs:
         with:
           asan: false
           ubsan: false
-          check_skip_long: false
 
   sanitizer_check:
     permissions:
@@ -46,7 +45,6 @@ jobs:
         with:
           asan: true
           ubsan: true
-          check_skip_long: true
       - name: Setting job link variable
         if: ${{ failure() }}
         run: |

doc/sphinx/installation.rst

@@ -108,20 +108,23 @@ CUDA SDK to make use of GPU computation:
 
     sudo apt install nvidia-cuda-toolkit
 
+Later in the installation instructions, you will see CMake commands of the
+form ``cmake ..`` with optional arguments, such as ``cmake .. -D WITH_CUDA=ON``
+to activate CUDA. These commands may need to be adapted depending on which
+operating system and CUDA version you are using.
+
 On Ubuntu 22.04, the default GCC compiler is too recent for nvcc and will fail
 to compile sources that rely on ``std::function``. You can either use GCC 10:
 
 .. code-block:: bash
 
     CC=gcc-10 CXX=g++-10 cmake .. -D WITH_CUDA=ON
-    make -j
 
 or alternatively install Clang 12 as a replacement for nvcc and GCC:
 
 .. code-block:: bash
 
     CC=clang-12 CXX=clang++-12 cmake .. -D WITH_CUDA=ON -D WITH_CUDA_COMPILER=clang
-    make -j
 
 On Ubuntu 20.04, the default GCC compiler is also too recent for nvcc and will
 generate compiler errors. You can either install an older version of GCC and

doc/sphinx/inter_non-bonded.rst

@@ -44,15 +44,13 @@ Tabulated interaction
 
     Feature ``TABULATED`` required.
 
-
 The interface for tabulated interactions are implemented in the
 :class:`~espressomd.interactions.TabulatedNonBonded` class. They can be configured
 via the following syntax::
 
   system.non_bonded_inter[type1, type2].tabulated.set_params(
       min='min', max='max', energy='energy', force='force')
 
-
 This defines an interaction between particles of the types ``type1`` and
 ``type2`` according to an arbitrary tabulated pair potential by linear interpolation.
 ``force`` specifies the tabulated forces and ``energy`` the energies as a function of the
@@ -646,19 +644,6 @@ and involved types.
 
 The samples folder contains the script :file:`/samples/drude_bmimpf6.py` with a
 fully polarizable, coarse grained ionic liquid where this approach is applied.
-To use the script, compile espresso with the following features:
-
-.. code-block:: c++
-
-    #define EXTERNAL_FORCES
-    #define MASS
-    #define THERMOSTAT_PER_PARTICLE
-    #define ROTATION
-    #define ROTATIONAL_INERTIA
-    #define ELECTROSTATICS
-    #define VIRTUAL_SITES_RELATIVE
-    #define LENNARD_JONES
-    #define THOLE
 
 .. _Anisotropic non-bonded interactions:
 

src/core/Particle.hpp

@@ -553,7 +553,6 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   }
   auto const &ext_force() const { return p.ext_force; }
   auto &ext_force() { return p.ext_force; }
-
 #else  // EXTERNAL_FORCES
   constexpr bool has_fixed_coordinates() const { return false; }
   constexpr bool is_fixed_along(int const) const { return false; }

src/core/constraints/Constraints.hpp

@@ -41,29 +41,29 @@ template <class ParticleRange, class Constraint> class Constraints {
 
 private:
   void reset_forces() const {
-    for (auto const &c : *this) {
-      c->reset_force();
+    for (auto const &constraint : *this) {
+      constraint->reset_force();
     }
   }
 
   container_type m_constraints;
 
 public:
-  void add(std::shared_ptr<Constraint> const &c) {
-    if (not c->fits_in_box(box_geo.length())) {
+  void add(std::shared_ptr<Constraint> const &constraint) {
+    if (not constraint->fits_in_box(box_geo.length())) {
       throw std::runtime_error("Constraint not compatible with box size.");
     }
-    assert(std::find(m_constraints.begin(), m_constraints.end(), c) ==
+    assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) ==
            m_constraints.end());
 
-    m_constraints.emplace_back(c);
+    m_constraints.emplace_back(constraint);
     on_constraint_change();
   }
-  void remove(std::shared_ptr<Constraint> const &c) {
-    assert(std::find(m_constraints.begin(), m_constraints.end(), c) !=
+  void remove(std::shared_ptr<Constraint> const &constraint) {
+    assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) !=
            m_constraints.end());
     m_constraints.erase(
-        std::remove(m_constraints.begin(), m_constraints.end(), c),
+        std::remove(m_constraints.begin(), m_constraints.end(), constraint),
         m_constraints.end());
     on_constraint_change();
   }
@@ -82,8 +82,8 @@ template <class ParticleRange, class Constraint> class Constraints {
     for (auto &p : particles) {
       auto const pos = folded_position(p.pos(), box_geo);
       ParticleForce force{};
-      for (auto const &c : *this) {
-        force += c->force(p, pos, t);
+      for (auto const &constraint : *this) {
+        force += constraint->force(p, pos, t);
       }
 
       p.f += force;

src/core/constraints/HomogeneousMagneticField.cpp

@@ -28,7 +28,7 @@ namespace Constraints {
 
 ParticleForce HomogeneousMagneticField::force(const Particle &p,
                                               const Utils::Vector3d &, double) {
-#if defined(ROTATION) && defined(DIPOLES)
+#ifdef DIPOLES
   return {{}, vector_product(p.calc_dip(), m_field)};
 #else
   return {};

src/core/grid_based_algorithms/lb_boundaries.cpp

@@ -52,8 +52,6 @@
 #include <stdexcept>
 #include <vector>
 
-using Utils::get_linear_index;
-
 namespace LBBoundaries {
 
 std::vector<std::shared_ptr<LBBoundary>> lbboundaries;
@@ -87,8 +85,8 @@ bool sanity_check_mach_limit() {
                      });
 }
 
-void ek_init_boundaries() {
-#if defined(CUDA) && defined(EK_BOUNDARIES)
+#if defined(EK_BOUNDARIES)
+static void ek_init_boundaries() {
   int number_of_boundnodes = 0;
 
   std::vector<float> host_wallcharge_species_density;
@@ -164,8 +162,8 @@ void ek_init_boundaries() {
     ek_init_species_density_wallcharge(host_wallcharge_species_density.data(),
                                        wallcharge_species);
   }
-#endif
 }
+#endif // defined(EK_BOUNDARIES)
 
 /** Initialize boundary conditions for all constraints in the system. */
 void lb_init_boundaries() {
@@ -175,7 +173,9 @@ void lb_init_boundaries() {
     }
 #if defined(CUDA)
 #if defined(LB_BOUNDARIES_GPU)
+#if defined(EK_BOUNDARIES)
     ek_init_boundaries();
+#endif
     unsigned number_of_boundnodes = 0;
     std::vector<int> host_boundary_node_list;
     std::vector<int> host_boundary_index_list;
@@ -237,6 +237,7 @@ void lb_init_boundaries() {
 #endif // defined (CUDA)
   } else if (lattice_switch == ActiveLB::CPU) {
 #if defined(LB_BOUNDARIES)
+    using Utils::get_linear_index;
     boost::for_each(lbfields, [](auto &f) { f.boundary = 0; });
 
     auto const node_pos = calc_node_pos(comm_cart);
@@ -285,7 +286,6 @@ REGISTER_CALLBACK(lb_collect_boundary_forces_local)
 
 Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) {
   Utils::Vector3d force{};
-#if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU)
   auto const it =
       boost::find_if(lbboundaries, [lbb](std::shared_ptr<LBBoundary> const &i) {
         return i.get() == lbb;
@@ -296,7 +296,7 @@ Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) {
                              "system.lbboundaries.");
   std::vector<double> forces(3 * lbboundaries.size());
   if (lattice_switch == ActiveLB::GPU) {
-#if defined(LB_BOUNDARIES_GPU) && defined(CUDA)
+#if defined(LB_BOUNDARIES_GPU)
     lb_gpu_get_boundary_forces(forces);
 #endif
   } else if (lattice_switch == ActiveLB::CPU) {
@@ -309,10 +309,9 @@ Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) {
   force[0] = forces[3 * container_index + 0];
   force[1] = forces[3 * container_index + 1];
   force[2] = forces[3 * container_index + 2];
-#endif
   return force;
 }
 
-#endif /* LB_BOUNDARIES or LB_BOUNDARIES_GPU */
+#endif // defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU)
 
 } // namespace LBBoundaries

src/core/magnetostatics/barnes_hut_gpu_cuda.cu

@@ -709,12 +709,12 @@ __global__ __launch_bounds__(THREADS5, FACTOR5) void forceCalculationKernel(
     // calculated within the array dq[i], which will
     // be compared later with squared distance between the particle
     // and the cell depending on a cell level.
-    // Original tree box edge (2*radiusd) should be divided *0.5
+    // Original tree box edge (2*radiusd) should be halved
     // as much as the tree depth takes place.
     dq[0] = radiusd * radiusd * *itolsqd;
     for (i = 1; i < maxdepthd; i++) {
-      dq[i] = dq[i - 1] * 0.25f;
-      dq[i - 1] += *epssqd;
+      dq[i] = dq[i - 1] * 0.25f; // halving of the squared distance
+      dq[i - 1] += *epssqd;      // increase thickness of previous cell
     }
     dq[i - 1] += *epssqd;
 
@@ -820,7 +820,7 @@ __global__ __launch_bounds__(THREADS5, FACTOR5) void forceCalculationKernel(
 #endif
               if (n != i) {
 
-                auto const d1 = sqrtf(tmp /*, 0.5f*/);
+                auto const d1 = sqrtf(tmp);
                 auto const dd5 = __fdividef(1.0f, tmp * tmp * d1);
                 auto b = 0.0f;
                 auto b2 = 0.0f;
@@ -891,7 +891,7 @@ __global__ __launch_bounds__(THREADS5, FACTOR5) void forceCalculationKernel(
 __global__ __launch_bounds__(THREADS5, FACTOR5) void energyCalculationKernel(
     float pf, float *energySum) {
   // NOTE: the algorithm of this kernel is almost identical to
-  // forceCalculationKernel. See comments there.
+  // @ref forceCalculationKernel. See comments there.
 
   int i, n, t;
   float dr[3], h[3], u[3], uc[3];
@@ -962,20 +962,16 @@ __global__ __launch_bounds__(THREADS5, FACTOR5) void energyCalculationKernel(
               dr[l] = -bhpara->r[3 * n + l] + bhpara->r[3 * i + l];
               tmp += dr[l] * dr[l];
             }
+            // check if all threads agree that cell is far enough away
+            // (or is a body)
 #if defined(__CUDACC__) && CUDA_VERSION >= 9000
             if ((n < bhpara->nbodies) ||
-                __all_sync(
-                    __activemask(),
-                    tmp >= dq[depth])) { // check if all threads agree that cell
-                                         // is far enough away (or is a body)
+                __all_sync(__activemask(), tmp >= dq[depth])) {
 #else
-            if ((n < bhpara->nbodies) ||
-                __all(tmp >=
-                      dq[depth])) { // check if all threads agree that cell
-                                    // is far enough away (or is a body)
+            if ((n < bhpara->nbodies) || __all(tmp >= dq[depth])) {
 #endif
               if (n != i) {
-                auto const d1 = sqrtf(tmp /*, 0.5f*/);
+                auto const d1 = sqrtf(tmp);
                 auto const dd5 = __fdividef(1.0f, tmp * tmp * d1);
                 auto b = 0.0f;
                 for (int l = 0; l < 3; l++) {
@@ -1275,4 +1271,4 @@ void fill_bh_data(float const *r, float const *dip, BHData const *bh_data) {
   cuda_safe_mem(cudaMemcpy(bh_data->u, dip, size, cudaMemcpyDeviceToDevice));
 }
 
-#endif // BARNES_HUT
+#endif // DIPOLAR_BARNES_HUT

src/core/magnetostatics/dds_gpu_cuda.cu

@@ -152,7 +152,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n,
   // There is one thread per particle. Each thread computes interactions
   // with particles whose id is smaller than the thread id.
   // The force and torque of all the interaction partners of the current thread
-  // is atomically added to global results ad once.
+  // is atomically added to global results at once.
   // The result for the particle id equal to the thread id is atomically added
   // to global memory at the end.
 

src/core/magnetostatics/dipoles_inline.hpp

@@ -51,7 +51,7 @@ struct ShortRangeForceKernel
     return {};
   }
 
-#ifdef P3M
+#ifdef DP3M
   result_type operator()(std::shared_ptr<DipolarP3M> const &ptr) const {
     auto const &actor = *ptr;
     return kernel_type{[&actor](Particle const &p1, Particle const &p2,
@@ -60,7 +60,7 @@ struct ShortRangeForceKernel
       return actor.pair_force(p1, p2, d, dist2, dist);
     }};
   }
-#endif // P3M
+#endif // DP3M
 
   result_type
   operator()(std::shared_ptr<DipolarLayerCorrection> const &ptr) const {
@@ -82,7 +82,7 @@ struct ShortRangeEnergyKernel
     return {};
   }
 
-#ifdef P3M
+#ifdef DP3M
   result_type operator()(std::shared_ptr<DipolarP3M> const &ptr) const {
     auto const &actor = *ptr;
     return kernel_type{[&actor](Particle const &p1, Particle const &p2,
@@ -91,7 +91,7 @@ struct ShortRangeEnergyKernel
       return actor.pair_energy(p1, p2, d, dist2, dist);
     }};
   }
-#endif // P3M
+#endif // DP3M
 
   result_type
   operator()(std::shared_ptr<DipolarLayerCorrection> const &ptr) const {

src/core/observables/ParticleTraits.hpp

@@ -41,7 +41,7 @@ template <> struct traits<Particle> {
   }
   auto charge(Particle const &p) const { return p.q(); }
   auto dipole_moment(Particle const &p) const {
-#if defined(ROTATION) && defined(DIPOLES)
+#ifdef DIPOLES
     return p.calc_dip();
 #else
     return Utils::Vector3d{};

src/core/p3m/common.hpp

@@ -39,17 +39,18 @@
 
 #include <utils/Vector.hpp>
 
+#include <array>
+#include <vector>
+
 /** This value indicates metallic boundary conditions. */
 auto constexpr P3M_EPSILON_METALLIC = 0.0;
 
 #if defined(P3M) || defined(DP3M)
 
 #include "LocalBox.hpp"
 
-#include <array>
 #include <cstddef>
 #include <stdexcept>
-#include <vector>
 
 namespace detail {
 /** @brief Index helpers for direct and reciprocal space.

src/core/p3m/fft.cpp

@@ -738,7 +738,7 @@ void fft_perform_back(double *data, bool check_complex, fft_data_struct &fft,
   for (int i = 0; i < fft.plan[1].new_size; i++) {
     fft.data_buf[i] = data[2 * i]; /* real value */
     // Vincent:
-    if (check_complex && (data[2 * i + 1] > 1e-5)) {
+    if (check_complex and std::abs(data[2 * i + 1]) > 1e-5) {
       printf("Complex value is not zero (i=%d,data=%g)!!!\n", i,
              data[2 * i + 1]);
       if (i > 100)