New variant implementation

CMakeLists.txt

@@ -805,11 +805,6 @@ endif()
 
 find_package(Boost 1.83.0 REQUIRED ${ESPRESSO_BOOST_COMPONENTS})
 
-# enable boost::variant with more than 20 types
-target_compile_options(
-  espresso_cpp_flags INTERFACE -DBOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
-                               -DBOOST_MPL_LIMIT_LIST_SIZE=30)
-
 #
 # Paths
 #

src/core/bond_breakage/actions.hpp

@@ -20,7 +20,6 @@
 #pragma once
 
 #include <boost/functional/hash.hpp>
-#include <boost/variant.hpp>
 
 #include <array>
 #include <cstddef>
@@ -74,7 +73,7 @@ struct DeleteAllBonds {
 } // namespace BondBreakage
 
 // Hash support for std::unordered_set
-namespace boost {
+namespace std {
 template <> struct hash<BondBreakage::DeleteBond> {
   std::size_t operator()(BondBreakage::DeleteBond const &t) const noexcept {
     return t.hash_value();
@@ -91,4 +90,4 @@ template <> struct hash<BondBreakage::DeleteAllBonds> {
     return t.hash_value();
   }
 };
-} // namespace boost
+} // namespace std

src/core/bond_breakage/bond_breakage.cpp

@@ -30,18 +30,18 @@
 
 #include <boost/mpi.hpp>
 #include <boost/serialization/access.hpp>
-#include <boost/variant.hpp>
 
 #include <cassert>
 #include <memory>
 #include <unordered_set>
 #include <utility>
+#include <variant>
 #include <vector>
 
 namespace BondBreakage {
 
 // Variant holding any of the actions
-using Action = boost::variant<DeleteBond, DeleteAngleBond, DeleteAllBonds>;
+using Action = std::variant<DeleteBond, DeleteAngleBond, DeleteAllBonds>;
 
 // Set of actions
 using ActionSet = std::unordered_set<Action>;
@@ -160,7 +160,7 @@ static void remove_pair_bonds_to(Particle &p, int other_pid) {
 }
 
 // Handler for the different delete events
-class execute : public boost::static_visitor<> {
+class execute {
   CellStructure &cell_structure;
 
 public:
@@ -199,7 +199,7 @@ void BondBreakage::process_queue_impl(System::System &system) {
 
   // Execute actions
   for (auto const &a : actions) {
-    boost::apply_visitor(execute(cell_structure), a);
+    std::visit(execute(cell_structure), a);
     system.on_particle_change();
   }
 }

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -27,31 +27,22 @@
 #include "thermalized_bond.hpp"
 #include "thermostat.hpp"
 
-#include <boost/variant.hpp>
-
 #include <algorithm>
 #include <numeric>
 #include <ranges>
-
-/** Visitor to get the bond cutoff from the bond parameter variant */
-class BondCutoff : public boost::static_visitor<double> {
-public:
-  template <typename T> double operator()(T const &bond) const {
-    return bond.cutoff();
-  }
-};
+#include <variant>
 
 double BondedInteractionsMap::maximal_cutoff() const {
   auto const max_cut_bonded = std::accumulate(
       begin(), end(), BONDED_INACTIVE_CUTOFF, [](auto max_cut, auto const &kv) {
-        return std::max(max_cut,
-                        boost::apply_visitor(BondCutoff(), *kv.second));
+        auto constexpr visitor = [](auto const &bond) { return bond.cutoff(); };
+        return std::max(max_cut, std::visit(visitor, *kv.second));
       });
 
   /* Check if there are dihedrals */
   auto const any_dihedrals = std::ranges::any_of(*this, [](auto const &kv) {
-    return (boost::get<DihedralBond>(&(*kv.second)) ||
-            boost::get<TabulatedDihedralBond>(&(*kv.second)));
+    return (std::holds_alternative<DihedralBond>(*kv.second) or
+            std::holds_alternative<TabulatedDihedralBond>(*kv.second));
   });
 
   /* dihedrals: the central particle is indirectly connected to the fourth
@@ -65,11 +56,11 @@ void BondedInteractionsMap::on_ia_change() {
   n_rigid_bonds = 0;
 #endif
   for (auto const &bond : std::views::elements<1>(*this)) {
-    if (boost::get<ThermalizedBond>(&(*bond)) != nullptr) {
+    if (std::holds_alternative<ThermalizedBond>(*bond)) {
       ++n_thermalized_bonds;
     }
 #ifdef BOND_CONSTRAINT
-    if (boost::get<RigidBond>(&(*bond)) != nullptr) {
+    if (std::holds_alternative<RigidBond>(*bond)) {
       ++n_rigid_bonds;
     }
 #endif
@@ -82,30 +73,30 @@ void BondedInteractionsMap::on_ia_change() {
 
 void BondedInteractionsMap::activate_bond(mapped_type const &ptr) {
   auto &system = get_system();
-  if (auto bond = boost::get<ThermalizedBond>(ptr.get())) {
+  if (auto bond = std::get_if<ThermalizedBond>(ptr.get())) {
     bond->set_thermostat_view(system.thermostat);
   }
-  if (auto bond = boost::get<IBMVolCons>(ptr.get())) {
+  if (auto bond = std::get_if<IBMVolCons>(ptr.get())) {
     system.immersed_boundaries->register_softID(*bond);
   }
-  if (auto bond = boost::get<IBMTriel>(ptr.get())) {
+  if (auto bond = std::get_if<IBMTriel>(ptr.get())) {
     bond->initialize(*system.box_geo, *system.cell_structure);
   }
-  if (auto bond = boost::get<IBMTribend>(ptr.get())) {
+  if (auto bond = std::get_if<IBMTribend>(ptr.get())) {
     bond->initialize(*system.box_geo, *system.cell_structure);
   }
 }
 
 void BondedInteractionsMap::deactivate_bond(mapped_type const &ptr) {
-  if (auto bond = boost::get<ThermalizedBond>(ptr.get())) {
+  if (auto bond = std::get_if<ThermalizedBond>(ptr.get())) {
     bond->unset_thermostat_view();
     n_thermalized_bonds = -1;
   }
-  if (auto bond = boost::get<IBMVolCons>(ptr.get())) {
+  if (auto bond = std::get_if<IBMVolCons>(ptr.get())) {
     bond->unset_volumes_view();
   }
 #ifdef BOND_CONSTRAINT
-  if (boost::get<RigidBond>(ptr.get()) != nullptr) {
+  if (std::get_if<RigidBond>(ptr.get())) {
     n_rigid_bonds = -1;
   }
 #endif

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -47,14 +47,13 @@
 #include "TabulatedPotential.hpp"
 #include "system/Leaf.hpp"
 
-#include <boost/variant.hpp>
-
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <optional>
 #include <stdexcept>
 #include <unordered_map>
+#include <variant>
 #include <vector>
 
 /* Special cutoff value for a disabled bond.
@@ -74,24 +73,16 @@ struct VirtualBond {
   double cutoff() const { return BONDED_INACTIVE_CUTOFF; }
 };
 
-/** Visitor to get the number of bound partners from the bond parameter
- *  variant.
- */
-class BondNumPartners : public boost::static_visitor<int> {
-public:
-  template <typename T> int operator()(T const &) const { return T::num; }
-};
-
 /** Variant in which to store the parameters of an individual bonded
  *  interaction
  */
 using Bonded_IA_Parameters =
-    boost::variant<NoneBond, FeneBond, HarmonicBond, QuarticBond, BondedCoulomb,
-                   BondedCoulombSR, AngleHarmonicBond, AngleCosineBond,
-                   AngleCossquareBond, DihedralBond, TabulatedDistanceBond,
-                   TabulatedAngleBond, TabulatedDihedralBond, ThermalizedBond,
-                   RigidBond, IBMTriel, IBMVolCons, IBMTribend,
-                   OifGlobalForcesBond, OifLocalForcesBond, VirtualBond>;
+    std::variant<NoneBond, FeneBond, HarmonicBond, QuarticBond, BondedCoulomb,
+                 BondedCoulombSR, AngleHarmonicBond, AngleCosineBond,
+                 AngleCossquareBond, DihedralBond, TabulatedDistanceBond,
+                 TabulatedAngleBond, TabulatedDihedralBond, ThermalizedBond,
+                 RigidBond, IBMTriel, IBMVolCons, IBMTribend,
+                 OifGlobalForcesBond, OifLocalForcesBond, VirtualBond>;
 
 /**
  * @brief container for bonded interactions.
@@ -147,7 +138,7 @@ class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
   auto size() const { return m_params.size(); }
   auto get_next_key() const { return next_key; }
   auto get_zero_based_type(int bond_id) const {
-    return contains(bond_id) ? at(bond_id)->which() : 0;
+    return contains(bond_id) ? static_cast<int>(at(bond_id)->index()) : 0;
   }
   auto get_n_thermalized_bonds() const {
     assert(n_thermalized_bonds >= 0);
@@ -198,7 +189,7 @@ class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
         bonds.begin(), bonds.end(),
         [this, partner_id = p_partner.id()](BondView const &bond) {
           auto const &bond_ptr = at(bond.bond_id());
-          return (boost::get<BondType>(bond_ptr.get()) != nullptr) and
+          return std::holds_alternative<BondType>(*bond_ptr.get()) and
                  (bond.partner_ids()[0] == partner_id);
         });
   }
@@ -234,5 +225,5 @@ class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
 
 /** @brief Get the number of bonded partners for the specified bond. */
 inline int number_of_partners(Bonded_IA_Parameters const &iaparams) {
-  return boost::apply_visitor(BondNumPartners(), iaparams);
+  return std::visit([]<typename T>(T const &) { return T::num; }, iaparams);
 }

src/core/bonded_interactions/bonded_interactions.dox

@@ -110,15 +110,6 @@
  *    - Include the header file containing the new bond type.
  *    - Add the new bond type to @ref Bonded_IA_Parameters at the end of the
  *      types list.
- *    - If by doing this, the length of the list in Bonded_IA_Parameters passes
- *      over a multiple of 10, you may have to update ESPResSo's top level
- *      CMakeLists.txt:
- *      @code
- *      # enable boost::variant with more than 20 types
- *      target_compile_options(
- *        espresso_cpp_flags INTERFACE -DBOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
- *                                     -DBOOST_MPL_LIMIT_LIST_SIZE=40)
- *      @endcode
  *  * In forces_inline.hpp:
  *    - A call to the new bond's force calculation needs to be placed in either
  *      of the functions @ref calc_bond_pair_force(), @ref
@@ -128,7 +119,7 @@
  *      example
  *      @code{.cpp}
  *      // ...
- *      else if (auto const *iap = boost::get<QuarticBond>(&iaparams)) {
+ *      else if (auto const *iap = std::get_if<QuarticBond>(&iaparams)) {
  *        return iap->force(dx);
  *      }
  *      // ...
@@ -141,7 +132,7 @@
  *      example
  *      @code{.cpp}
  *      // ...
- *      else if (auto const *iap = boost::get<QuarticBond>(&iaparams)) {
+ *      else if (auto const *iap = std::get_if<QuarticBond>(&iaparams)) {
  *        return iap->energy(dx);
  *      }
  *      // ...

src/core/cell_system/CellStructure.cpp

@@ -40,7 +40,6 @@
 #include <utils/math/sqr.hpp>
 
 #include <boost/mpi/collectives/all_reduce.hpp>
-#include <boost/variant.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -53,6 +52,7 @@
 #include <stdexcept>
 #include <string>
 #include <utility>
+#include <variant>
 #include <vector>
 
 #ifdef SHARED_MEMORY_PARALLELISM
@@ -231,7 +231,7 @@ void CellStructure::resort_particles(bool global_flag) {
   m_decomposition->resort(global_flag, diff);
 
   for (auto d : diff) {
-    boost::apply_visitor(UpdateParticleIndexVisitor{this}, d);
+    std::visit(UpdateParticleIndexVisitor{this}, d);
   }
 
   auto const &lebc = get_system().box_geo->lees_edwards_bc();

src/core/cell_system/ParticleDecomposition.hpp

@@ -26,10 +26,9 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/variant.hpp>
-
 #include <optional>
 #include <span>
+#include <variant>
 #include <vector>
 
 struct RemovedParticle {
@@ -43,7 +42,7 @@ struct ModifiedList {
 /**
  * @brief Change of Particle Address.
  */
-using ParticleChange = boost::variant<RemovedParticle, ModifiedList>;
+using ParticleChange = std::variant<RemovedParticle, ModifiedList>;
 
 /**
  * @brief A distributed particle decomposition.

src/core/communication.hpp

@@ -28,9 +28,9 @@
  *  The asynchronous MPI communication is used during the script
  *  evaluation. Except for the head node that interprets the interface
  *  script, all other nodes wait in @ref mpi_loop() for the head node to
- *  issue an action using @ref mpi_call(). @ref mpi_loop() immediately
+ *  issue an action using @c MpiCallbacks::call(). @ref mpi_loop() immediately
  *  executes an @c MPI_Bcast and therefore waits for the head node to
- *  broadcast a command, which is done by @ref mpi_call(). The request
+ *  broadcast a command, which is done by @c MpiCallbacks::call(). The request
  *  consists of a callback function with an arbitrary number of arguments.
  *
  *  To add new actions (e.g. to implement new interface functionality), do the