Add support for ARM architectures (#4708)

Description of changes:
- adjust test tolerances for architectures that don't have double extended precision, e.g. ARM
- add sanity checks in chain analysis functions and handle the corner case `number_of_chains=1`
- re-enable support for fast-math mode

maintainer/CI/build_cmake.sh

@@ -141,7 +141,7 @@ if [ "${with_coverage}" = true ]; then
 fi
 
 if [ "${with_fast_math}" = true ]; then
-    cmake_param_protected="-DCMAKE_CXX_FLAGS=-ffast-math -fno-finite-math-only"
+    cmake_param_protected="-DCMAKE_CXX_FLAGS=-ffast-math"
 fi
 
 cmake_params="-D CMAKE_BUILD_TYPE=${build_type} -D CMAKE_CXX_STANDARD=${with_cxx_standard} -D ESPRESSO_WARNINGS_ARE_ERRORS=ON ${cmake_params}"

src/core/analysis/statistics_chain.cpp

@@ -35,12 +35,11 @@
 #include <cmath>
 #include <stdexcept>
 
-std::array<double, 4> calc_re(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 4> calc_re(int chain_start, int n_chains, int chain_length) {
   double dist = 0.0, dist2 = 0.0, dist4 = 0.0;
   std::array<double, 4> re;
 
-  for (int i = 0; i < chain_n_chains; i++) {
+  for (int i = 0; i < n_chains; i++) {
     auto const &p1 =
         get_particle_data(chain_start + i * chain_length + chain_length - 1);
     auto const &p2 = get_particle_data(chain_start + i * chain_length);
@@ -53,20 +52,19 @@ std::array<double, 4> calc_re(int chain_start, int chain_n_chains,
     dist2 += norm2;
     dist4 += norm2 * norm2;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(n_chains);
   re[0] = dist / tmp;
   re[2] = dist2 / tmp;
-  re[1] = sqrt(re[2] - re[0] * re[0]);
-  re[3] = sqrt(dist4 / tmp - re[2] * re[2]);
+  re[1] = (n_chains == 1) ? 0. : std::sqrt(re[2] - Utils::sqr(re[0]));
+  re[3] = (n_chains == 1) ? 0. : std::sqrt(dist4 / tmp - Utils::sqr(re[2]));
   return re;
 }
 
-std::array<double, 4> calc_rg(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 4> calc_rg(int chain_start, int n_chains, int chain_length) {
   double r_G = 0.0, r_G2 = 0.0, r_G4 = 0.0;
   std::array<double, 4> rg;
 
-  for (int i = 0; i < chain_n_chains; i++) {
+  for (int i = 0; i < n_chains; i++) {
     double M = 0.0;
     Utils::Vector3d r_CM{};
     for (int j = 0; j < chain_length; j++) {
@@ -94,22 +92,21 @@ std::array<double, 4> calc_rg(int chain_start, int chain_n_chains,
     r_G2 += tmp;
     r_G4 += tmp * tmp;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(n_chains);
   rg[0] = r_G / tmp;
   rg[2] = r_G2 / tmp;
-  rg[1] = sqrt(rg[2] - Utils::sqr(rg[0]));
-  rg[3] = sqrt(r_G4 / tmp - Utils::sqr(rg[2]));
+  rg[1] = (n_chains == 1) ? 0. : std::sqrt(rg[2] - Utils::sqr(rg[0]));
+  rg[3] = (n_chains == 1) ? 0. : std::sqrt(r_G4 / tmp - Utils::sqr(rg[2]));
   return rg;
 }
 
-std::array<double, 2> calc_rh(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 2> calc_rh(int chain_start, int n_chains, int chain_length) {
   double r_H = 0.0, r_H2 = 0.0;
   std::array<double, 2> rh;
 
   auto const chain_l = static_cast<double>(chain_length);
   auto const prefac = 0.5 * chain_l * (chain_l - 1.);
-  for (int p = 0; p < chain_n_chains; p++) {
+  for (int p = 0; p < n_chains; p++) {
     double ri = 0.0;
     for (int i = chain_start + chain_length * p;
          i < chain_start + chain_length * (p + 1); i++) {
@@ -126,8 +123,8 @@ std::array<double, 2> calc_rh(int chain_start, int chain_n_chains,
     r_H += tmp;
     r_H2 += tmp * tmp;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(n_chains);
   rh[0] = r_H / tmp;
-  rh[1] = sqrt(r_H2 / tmp - Utils::sqr(rh[0]));
+  rh[1] = (n_chains == 1) ? 0. : std::sqrt(r_H2 / tmp - Utils::sqr(rh[0]));
   return rh;
 }

src/core/analysis/statistics_chain.hpp

@@ -34,11 +34,10 @@
  * of monodisperse polymers with continuous ids.
  *
  * @param chain_start The id of the first monomer of the first chain.
- * @param chain_n_chains Number of chains contained in the range.
+ * @param n_chains Number of chains contained in the range.
  * @param chain_length The length of every chain.
  */
-std::array<double, 4> calc_re(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 4> calc_re(int chain_start, int n_chains, int chain_length);
 
 /**
  * @brief Calculate the radius of gyration.
@@ -47,11 +46,10 @@ std::array<double, 4> calc_re(int chain_start, int chain_n_chains,
  * of monodisperse polymers with continuous ids.
  *
  * @param chain_start The id of the first monomer of the first chain.
- * @param chain_n_chains Number of chains contained in the range.
+ * @param n_chains Number of chains contained in the range.
  * @param chain_length The length of every chain.
  */
-std::array<double, 4> calc_rg(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 4> calc_rg(int chain_start, int n_chains, int chain_length);
 
 /**
  * @brief Calculate the hydrodynamic radius (ref. Kirkwood-Zimm theory).
@@ -60,10 +58,9 @@ std::array<double, 4> calc_rg(int chain_start, int chain_n_chains,
  * of monodisperse polymers with continuous ids.
  *
  * @param chain_start The id of the first monomer of the first chain.
- * @param chain_n_chains Number of chains contained in the range.
+ * @param n_chains Number of chains contained in the range.
  * @param chain_length The length of every chain.
  */
-std::array<double, 2> calc_rh(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 2> calc_rh(int chain_start, int n_chains, int chain_length);
 
 #endif

src/core/electrostatics/icc.cpp

@@ -129,6 +129,12 @@ void ICCStar::iteration(CellStructure &cell_structure,
     for (auto &p : particles) {
       auto const pid = p.id();
       if (pid >= icc_cfg.first_id and pid < icc_cfg.n_icc + icc_cfg.first_id) {
+        if (p.q() == 0.) {
+          runtimeErrorMsg()
+              << "ICC found zero electric charge on a particle. This must "
+                 "never happen";
+          break;
+        }
         auto const id = p.id() - icc_cfg.first_id;
         /* the dielectric-related prefactor: */
         auto const eps_in = icc_cfg.epsilons[id];
@@ -176,7 +182,7 @@ void ICCStar::iteration(CellStructure &cell_structure,
               << "Particle with id " << p.id() << " has a charge (q=" << p.q()
               << ") that is too large for the ICC algorithm";
 
-          max_rel_diff = std::numeric_limits<double>::infinity();
+          max_rel_diff = std::numeric_limits<double>::max();
           break;
         }
       }

src/core/particle_node.cpp

@@ -457,18 +457,6 @@ static bool maybe_move_particle(int p_id, Utils::Vector3d const &pos) {
   return true;
 }
 
-void particle_checks(int p_id, Utils::Vector3d const &pos) {
-  if (p_id < 0) {
-    throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
-  }
-#ifndef __FAST_MATH__
-  if (std::isnan(pos[0]) or std::isnan(pos[1]) or std::isnan(pos[2]) or
-      std::isinf(pos[0]) or std::isinf(pos[1]) or std::isinf(pos[2])) {
-    throw std::domain_error("Particle position must be finite");
-  }
-#endif // __FAST_MATH__
-}
-
 void remove_all_particles() {
   ::cell_structure.remove_all_particles();
   on_particle_change();

src/core/particle_node.hpp

@@ -82,7 +82,6 @@ void clear_particle_node();
  * @param pos   The particle position.
  */
 void make_new_particle(int p_id, Utils::Vector3d const &pos);
-void particle_checks(int p_id, Utils::Vector3d const &pos);
 
 /**
  * @brief Move particle to a new position.

src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp

@@ -68,7 +68,7 @@ class ReactionAlgorithm : public ReactionMethods::ReactionAlgorithm {
 // Check the base class for all Monte Carlo algorithms.
 BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
   using ReactionMethods::SingleReaction;
-  auto constexpr tol = 100. * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   auto const comm = boost::mpi::communicator();
 
   // check acceptance rate

src/core/reaction_methods/tests/SingleReaction_test.cpp

@@ -32,7 +32,7 @@
 // and the configurational move probability for a given system state.
 BOOST_AUTO_TEST_CASE(SingleReaction_test) {
   using namespace ReactionMethods;
-  constexpr double tol = 100 * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
   // create a reaction A -> 3 B + 4 C
   int const type_A = 0;
@@ -64,7 +64,7 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) {
             std::unordered_map<int, int>{{type_A, i}, {type_B, j}, {type_C, k}};
         auto const val = calculate_factorial_expression(reaction, p_numbers);
         auto const ref = g(i, -1) * g(j, 3) * g(k, 4);
-        BOOST_CHECK_CLOSE(val, ref, 5 * tol);
+        BOOST_CHECK_CLOSE(val, ref, 5. * tol);
       }
     }
   }

src/core/reaction_methods/tests/reaction_methods_utils_test.cpp

@@ -28,7 +28,7 @@
 
 BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) {
   using namespace ReactionMethods;
-  constexpr double tol = 100. * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
   auto const reaction_ensemble_combinations = [](int N, int nu) {
     return (N + nu < 0) ? 0. : std::tgamma(N + 1) / std::tgamma(N + nu + 1);

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -92,7 +92,7 @@ static auto copy_particle_to_head_node(boost::mpi::communicator const &comm,
 }
 
 BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
-  constexpr auto tol = 100. * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   auto const comm = boost::mpi::communicator();
   auto const rank = comm.rank();
   auto const n_nodes = comm.size();

src/core/unit_tests/Verlet_list_test.cpp

@@ -167,7 +167,7 @@ auto const propagators =
 BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update,
                        bdata::make(node_grids) * bdata::make(propagators),
                        node_grid, integration_helper) {
-  constexpr auto tol = 100. * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   auto const comm = boost::mpi::communicator();
   auto const rank = comm.rank();
 

src/core/unit_tests/field_coupling_couplings_test.cpp

@@ -16,7 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#define BOOST_TEST_MODULE AutoParameter test
+
+#define BOOST_TEST_MODULE Field coupling test for fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 

src/core/unit_tests/field_coupling_fields_test.cpp

@@ -17,7 +17,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#define BOOST_TEST_MODULE AutoParameter test
+#define BOOST_TEST_MODULE Field coupling test for fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 
@@ -37,7 +37,7 @@
 #include <limits>
 #include <type_traits>
 
-constexpr auto eps = 10 * std::numeric_limits<double>::epsilon();
+auto constexpr eps = 10. * std::numeric_limits<double>::epsilon();
 
 using namespace FieldCoupling::Fields;
 
@@ -269,6 +269,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) {
     BOOST_CHECK_SMALL(std::abs(interpolated_value - field_value), eps);
   }
 
+#ifndef __FAST_MATH__
   /* jacobian value */
   {
     using Utils::Interpolation::bspline_3d_gradient_accumulate;
@@ -296,6 +297,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) {
 
     BOOST_CHECK_SMALL((interpolated_value - field_value).norm(), eps);
   }
+#endif // __FAST_MATH__
 }
 
 BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
@@ -342,6 +344,7 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
     BOOST_CHECK_SMALL((interpolated_value - field_value).norm(), eps);
   }
 
+#ifndef __FAST_MATH__
   /* jacobian value */
   {
     using Utils::Interpolation::bspline_3d_gradient_accumulate;
@@ -380,4 +383,5 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) {
     BOOST_CHECK_SMALL(
         (interpolated_value.row<1>() - field_value.row<1>()).norm(), eps);
   }
+#endif // __FAST_MATH__
 }

src/core/unit_tests/field_coupling_force_field_test.cpp

@@ -16,7 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#define BOOST_TEST_MODULE test
+
+#define BOOST_TEST_MODULE Field coupling test for force fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 

src/core/unit_tests/p3m_test.cpp

@@ -64,7 +64,7 @@ BOOST_AUTO_TEST_CASE(calc_meshift_true) {
 #if defined(P3M) || defined(DP3M)
 BOOST_AUTO_TEST_CASE(analytic_cotangent_sum) {
   auto constexpr kernel = p3m_analytic_cotangent_sum;
-  auto constexpr tol = 100. * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
   // check only trivial cases
   for (auto const cao : {1, 2, 3, 4, 5, 6, 7}) {

src/core/unit_tests/periodic_fold_test.cpp

@@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(with_image_count) {
     BOOST_CHECK(res.first >= 0.);
     BOOST_CHECK(res.first <= box);
     BOOST_CHECK(std::abs(res.first - x + res.second * box) <=
-                std::numeric_limits<double>::epsilon());
+                4. * std::numeric_limits<double>::epsilon());
   }
 
   /* Corner right */
@@ -117,7 +117,7 @@ BOOST_AUTO_TEST_CASE(with_image_count) {
     BOOST_CHECK(res.first >= 0.);
     BOOST_CHECK(res.first < box);
     BOOST_CHECK(std::abs(res.first - x + res.second * box) <=
-                std::numeric_limits<double>::epsilon());
+                4. * std::numeric_limits<double>::epsilon());
   }
 }
 

src/core/unit_tests/specfunc_test.cpp

@@ -26,7 +26,7 @@
 #include <cmath>
 #include <limits>
 
-constexpr auto eps = 100. * std::numeric_limits<double>::epsilon();
+auto constexpr eps = 8. * 100. * std::numeric_limits<double>::epsilon();
 
 BOOST_AUTO_TEST_CASE(hurwitz_zeta_function) {
   constexpr auto max_bits = 54.0;

src/core/unit_tests/thermostats_test.cpp

@@ -41,8 +41,8 @@
 #include <tuple>
 
 // multiply by 100 because BOOST_CHECK_CLOSE takes a percentage tolerance,
-// and by 6 to account for error accumulation in thermostat functions
-constexpr auto tol = 6 * 100 * std::numeric_limits<double>::epsilon();
+// and by 8 to account for error accumulation in thermostat functions
+auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
 Particle particle_factory() {
   Particle p{};

src/particle_observables/tests/algorithms.cpp

@@ -79,7 +79,7 @@ BOOST_AUTO_TEST_CASE(algorithms_integer) {
 }
 
 BOOST_AUTO_TEST_CASE(algorithms_double) {
-  auto constexpr tol = 100 * std::numeric_limits<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   std::vector<double> const values{1., 2., 3., 4.};
   {
     auto const res = WeightedAverage<Testing::Identity, Testing::One>()(values);

src/script_interface/CMakeLists.txt

@@ -60,6 +60,10 @@ target_link_libraries(
   PUBLIC espresso::utils MPI::MPI_CXX Boost::mpi espresso::shapes
   PRIVATE espresso::cpp_flags)
 
+set_source_files_properties(
+  ${CMAKE_CURRENT_SOURCE_DIR}/particle_data/ParticleHandle.cpp
+  PROPERTIES COMPILE_FLAGS -fno-finite-math-only)
+
 target_include_directories(espresso_script_interface
                            PUBLIC ${CMAKE_SOURCE_DIR}/src)
 

src/script_interface/analysis/Analysis.cpp

@@ -48,6 +48,12 @@ namespace Analysis {
 
 /** @brief Check if a contiguous range of particle ids exists. */
 static void check_topology(int chain_start, int chain_length, int n_chains) {
+  if (n_chains <= 0) {
+    throw std::domain_error("Chain analysis needs at least 1 chain");
+  }
+  if (chain_length <= 0) {
+    throw std::domain_error("Chain analysis needs at least 1 bead per chain");
+  }
   for (int i = 0; i < chain_length * n_chains; ++i) {
     auto const pid = chain_start + i;
     if (not particle_exists(pid)) {