Fix issues with cluster analysis and chain analysis

doc/sphinx/analysis.rst

@@ -266,65 +266,45 @@ Chains
 ~~~~~~
 
 All analysis functions in this section require the topology of the chains to be set correctly.
-The above set of functions is designed to facilitate analysis of molecules.
-Molecules are expected to be a group of particles comprising a contiguous range of particle IDs.
-Each molecule is a set of consecutively numbered particles and all molecules are supposed to consist of the same number of particles.
-
-Some functions in this group require that the particles constituting a molecule are connected into
-linear chains (particle :math:`n` is connected to :math:`n+1` and so on)
-while others are applicable to molecules of whatever topology.
-
-
-.. _End to end distance:
-
-End-to-end distance
-^^^^^^^^^^^^^^^^^^^
-:meth:`espressomd.analyze.Analysis.calc_re`
-
-Returns the quadratic end-to-end-distance and its root averaged over all chains.
-
-.. _Radius of gyration:
-
-Radius of gyration
-^^^^^^^^^^^^^^^^^^
-:meth:`espressomd.analyze.Analysis.calc_rg`
-
-Returns the radius of gyration averaged over all chains.
-It is a radius of a sphere, which would have the same moment of inertia as the
-molecule, defined as
-
-.. math::
-
-   \label{eq:Rg}
-   R_{\mathrm G}^2 = \frac{1}{N} \sum\limits_{i=1}^{N} \left(\vec r_i - \vec r_{\mathrm{cm}}\right)^2\,,
-
-where :math:`\vec r_i` are position vectors of individual particles
-constituting a molecule and :math:`\vec r_{\mathrm{cm}}` is the position
-vector of its center of mass. The sum runs over all :math:`N` particles
-comprising the molecule. For more information see any polymer science
-book, e.g. :cite:`rubinstein03a`.
-
-
-.. _Hydrodynamic radius:
-
-Hydrodynamic radius
-^^^^^^^^^^^^^^^^^^^
-:meth:`espressomd.analyze.Analysis.calc_rh`
-
-Returns the hydrodynamic radius averaged over all chains.
-The following formula is used for the computation:
-
-.. math::
-
-   \label{eq:Rh}
-   \frac{1}{R_{\mathrm H}} = \frac{2}{N(N-1)} \sum\limits_{i=1}^{N} \sum\limits_{j<i}^{N} \frac{1}{|\vec r_i - \vec r_j|}\,,
-
-The above-mentioned formula is only valid under certain assumptions.
-For more information, see chapter 4 and equation 4.102 in :cite:`doi86a`.
-Note that the hydrodynamic radius is sometimes defined in a similar fashion
-but with a denominator of :math:`N^2` instead of :math:`N(N-1)` in the prefactor.
-Both versions are equivalent in the :math:`N\rightarrow \infty` limit but give
-numerically different values for finite polymers.
+A chain needs to be set up with a contiguous range of particle IDs, and the
+head resp. tail particle must be have the first resp. last id in the range.
+For chain observables that rely on connectivity information, the chain topology
+must be linear, i.e. particle :math:`n` is connected to :math:`n+1` and so on.
+Each chain is a set of consecutively numbered particles and all chains are
+supposed to consist of the same number of particles.
+
+The particles :attr:`~espressomd.particle_data.ParticleHandle.image_box`
+values must also be consistent, since these observables are calculated using
+*unfolded coordinates*. For example, if all particles were to be inserted in
+the central box except for the tail particle, which would be e.g. inserted in
+the fourth periodic image, the end-to-end distance and radius of gyration
+would be quite large, even if the tail particle is in close proximity to the
+penultimate particle in *folded coordinates*, because the last inter-particle
+distance would be evaluated as being 4 times the box length plus the bond
+length. Particles *can* have different ``image_box`` values, for example
+in a chain with equilibrium bond length 2, if the penultimate particle is at
+``[box_l - 1, 0, 0]`` and the tail particle at ``[box_l + 1, 0, 0]``, their
+inter-particle distance would be evaluated as 2 and the observables would
+be correctly calculated. This can lead to counter-intuitive behavior when
+chains are longer than half the box size in a fully periodic system. As an
+example, consider the end-to-end distance of a linear polymer growing on
+the x-axis. While :meth:`espressomd.system.System.distance()` would feature
+a triangle wave with a period equal to the box length in the x-direction,
+:meth:`espressomd.analyze.Analysis.calc_re` would be monotonically increasing,
+assuming the ``image_box`` values were properly set up. This is important to
+consider when setting up systems where the polymers are much longer than the
+box length, which is common when simulating polymer melts.
+
+.. _Available chain analysis functions:
+
+Available chain analysis functions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+* :meth:`espressomd.analyze.Analysis.calc_re`: average end-to-end-distance
+
+* :meth:`espressomd.analyze.Analysis.calc_rg`: average radius of gyration
+
+* :meth:`espressomd.analyze.Analysis.calc_rh`: average hydrodynamic radius
 
 
 .. _Observables framework:

src/core/BoxGeometry.hpp

@@ -188,12 +188,6 @@ class BoxGeometry {
                                     const Utils::Vector<T, 3> &b) const {
     if (type() == BoxType::LEES_EDWARDS) {
       auto const &shear_plane_normal = lees_edwards_bc().shear_plane_normal;
-      //      if (a[shear_plane_normal] <0 or a[shear_plane_normal]
-      //      >=length()[shear_plane_normal]
-      //          or b[shear_plane_normal] <0 or b[shear_plane_normal]
-      //          >=length()[shear_plane_normal]) throw
-      //          std::runtime_error("Unfolded position in shear plane normal
-      //          direction in LE get_mi_vector will lead to wrong result");
       auto a_tmp = a;
       auto b_tmp = b;
       a_tmp[shear_plane_normal] = Algorithm::periodic_fold(

src/core/bonded_interactions/bonded_interactions.dox

@@ -60,10 +60,12 @@
  *    double cutoff() const { return r0 + drmax; }
  *    @endcode
  *    The return value of \c cutoff() should be as large as the interaction
- *    range of the new interaction. This is only relevant to pairwise bonds
- *    and dihedral bonds. In all other cases, the return value should be -1,
- *    namely angle bonds as well as other bonds that don't have an interaction
- *    range.
+ *    range of the new interaction. This is only relevant to pairwise bonds.
+ *    In all other cases, the return value should be 0, namely angle bonds,
+ *    dihedral bonds as well as other bonds that don't have an interaction
+ *    range. The @ref VirtualBond is the exception to this rule; its range
+ *    is @ref BONDED_INACTIVE_CUTOFF to ensure that it is always skipped by
+ *    the short-range loop.
  *  * @code{.cpp}
  *    boost::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;
  *    @endcode

src/core/cluster_analysis/Cluster.cpp

@@ -37,6 +37,7 @@
 #include <cmath>
 #include <cstddef>
 #include <numeric>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -50,6 +51,7 @@ Utils::Vector3d Cluster::center_of_mass() {
 // Center of mass of an aggregate
 Utils::Vector3d
 Cluster::center_of_mass_subcluster(std::vector<int> const &particle_ids) {
+  sanity_checks();
   Utils::Vector3d com{};
 
   // The distances between the particles are "folded", such that all distances
@@ -79,6 +81,7 @@ Cluster::center_of_mass_subcluster(std::vector<int> const &particle_ids) {
 }
 
 double Cluster::longest_distance() {
+  sanity_checks();
   double ld = 0.;
   for (auto a = particles.begin(); a != particles.end(); a++) {
     for (auto b = a; ++b != particles.end();) {
@@ -101,6 +104,7 @@ double Cluster::radius_of_gyration() {
 
 double
 Cluster::radius_of_gyration_subcluster(std::vector<int> const &particle_ids) {
+  sanity_checks();
   // Center of mass
   Utils::Vector3d com = center_of_mass_subcluster(particle_ids);
   double sum_sq_dist = 0.;
@@ -128,6 +132,7 @@ std::vector<std::size_t> sort_indices(const std::vector<T> &v) {
 
 std::pair<double, double> Cluster::fractal_dimension(double dr) {
 #ifdef GSL
+  sanity_checks();
   Utils::Vector3d com = center_of_mass();
   // calculate Df using linear regression on the logarithms of the radii of
   // gyration against the number of particles in sub-clusters. Particles are
@@ -176,4 +181,11 @@ std::pair<double, double> Cluster::fractal_dimension(double dr) {
 #endif
 }
 
+void Cluster::sanity_checks() const {
+  if (::box_geo.type() != BoxType::CUBOID) {
+    throw std::runtime_error(
+        "Cluster analysis is not compatible with non-cuboid box types");
+  }
+}
+
 } // namespace ClusterAnalysis

src/core/cluster_analysis/Cluster.hpp

@@ -53,6 +53,9 @@ class Cluster {
    *
    *  @return fractal dimension, rms error of the fit */
   std::pair<double, double> fractal_dimension(double dr);
+
+private:
+  void sanity_checks() const;
 };
 
 } // namespace ClusterAnalysis

src/core/cluster_analysis/ClusterStructure.cpp

@@ -18,16 +18,19 @@
  */
 #include "ClusterStructure.hpp"
 
+#include "BoxGeometry.hpp"
 #include "Cluster.hpp"
 #include "PartCfg.hpp"
 #include "errorhandling.hpp"
+#include "grid.hpp"
 #include "partCfg_global.hpp"
 #include "particle_node.hpp"
 
 #include <utils/for_each_pair.hpp>
 
 #include <algorithm>
 #include <memory>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -49,6 +52,7 @@ inline bool ClusterStructure::part_of_cluster(const Particle &p) {
 void ClusterStructure::run_for_all_pairs() {
   // clear data structs
   clear();
+  sanity_checks();
 
   // Iterate over pairs
   Utils::for_each_pair(partCfg().begin(), partCfg().end(),
@@ -60,6 +64,7 @@ void ClusterStructure::run_for_all_pairs() {
 
 void ClusterStructure::run_for_bonded_particles() {
   clear();
+  sanity_checks();
   for (const auto &p : partCfg()) {
     for (auto const bond : p.bonds()) {
       if (bond.partner_ids().size() == 1) {
@@ -189,4 +194,11 @@ int ClusterStructure::get_next_free_cluster_id() {
   return max_seen_cluster + 1;
 }
 
+void ClusterStructure::sanity_checks() const {
+  if (::box_geo.type() != BoxType::CUBOID) {
+    throw std::runtime_error(
+        "Cluster analysis is not compatible with non-cuboid box types");
+  }
+}
+
 } // namespace ClusterAnalysis

src/core/cluster_analysis/ClusterStructure.hpp

@@ -76,6 +76,7 @@ class ClusterStructure {
   inline int find_id_for(int x);
   /** @brief Get next free cluster id */
   inline int get_next_free_cluster_id();
+  void sanity_checks() const;
 };
 
 } // namespace ClusterAnalysis

src/core/virtual_sites/VirtualSitesRelative.cpp

@@ -116,17 +116,6 @@ void VirtualSitesRelative::update() const {
 
     auto const &p_ref = *p_ref_ptr;
     p.pos() = p_ref.pos() + connection_vector(p_ref, p);
-    //    auto new_pos = p_ref.pos() + connection_vector(p_ref, p);
-    //    /* The shift has to respect periodic boundaries: if the reference
-    //     * particles is not in the same image box, we potentially avoid
-    //     shifting
-    //     * to the other side of the box. */
-    //    auto shift = box_geo.get_mi_vector(new_pos, p.pos());
-    //    p.pos() += shift;
-    //    Utils::Vector3i image_shift{};
-    //    fold_position(shift, image_shift, box_geo);
-    //    p.image_box() = p_ref.image_box() - image_shift;
-    //
     p.v() = velocity(p_ref, p);
 
     if (box_geo.type() == BoxType::LEES_EDWARDS) {

src/python/espressomd/analyze.py

@@ -237,6 +237,19 @@ class Analysis(ScriptInterfaceHelper):
         numbered, the last particle in that topology having id number
         ``chain_start + number_of_chains * chain_length - 1``.
 
+        The radius of gyration is the radius of a sphere which would have
+        the same moment of inertia as a polymer, and is defined as
+
+        .. math::
+
+           R_{\\mathrm G}^2 = \\frac{1}{N} \\sum\\limits_{i=1}^{N} \\left(\\vec r_i - \\vec r_{\\mathrm{cm}}\\right)^2\\,,
+
+        where :math:`\\vec r_i` are position vectors of individual particles
+        constituting the polymer and :math:`\\vec r_{\\mathrm{cm}}` is the
+        position of its center of mass. The sum runs over all :math:`N`
+        particles comprising the polymer. For more information see any
+        polymer science book, e.g. :cite:`rubinstein03a`.
+
         Parameters
         ----------
         chain_start : :obj:`int`
@@ -254,14 +267,28 @@ class Analysis(ScriptInterfaceHelper):
             its standard deviation.
 
     calc_rh()
-        Calculate the hydrodynamic mean radius of chains and its standard
+        Calculate the mean hydrodynamic radius of chains and its standard
         deviation.
 
         This requires that a set of chains of equal length which start
         with the particle number ``chain_start`` and are consecutively
         numbered, the last particle in that topology having id number
         ``chain_start + number_of_chains * chain_length - 1``.
 
+        The following formula is used for the calculation:
+
+        .. math::
+
+           \\frac{1}{R_{\\mathrm H}} = \\frac{2}{N(N-1)} \\sum\\limits_{i=1}^{N} \\sum\\limits_{j<i}^{N} \\frac{1}{|\\vec r_i - \\vec r_j|}\\,,
+
+        This formula is only valid under certain assumptions. For more
+        information, see chapter 4 and equation 4.102 in :cite:`doi86a`.
+        Note that the hydrodynamic radius is sometimes defined in a similar
+        fashion but with a denominator of :math:`N^2` instead of :math:`N(N-1)`
+        in the prefactor. Both versions are equivalent in the
+        :math:`N\\rightarrow \\infty` limit but give numerically different
+        values for finite polymers.
+
         Parameters
         ----------
         chain_start : :obj:`int`