Cluster analysis methods don't work with Lees-Edwards

[jngrad]

The particle position folding code in the cluster analysis methods don't support Lees-Edwards boundary conditions.

Below is a MWE created from test cases analyze_chains.py and cluster_analysis.py. The test case sets a polymer in such a way that all particles are on a straight line with a constant pairwise distance. All cluster methods seem to apply the shear offset twice when unfolding the coordinates. The GSL method doesn't properly unfold the coordinates at all.

To run this MWE, you first need to do git revert 47cb7a586.

import unittest as ut
import numpy as np
import espressomd
import espressomd.interactions
import espressomd.lees_edwards
import espressomd.pair_criteria
import espressomd.cluster_analysis

class Test(ut.TestCase):

    system = espressomd.System(box_l=(1, 1, 1))
    system.periodicity = [True, True, True]
    system.time_step = 0.01
    system.cell_system.skin = 0.1

    fene = espressomd.interactions.FeneBond(k=1, d_r_max=0.05)
    system.bonded_inter.add(fene)

    cs = espressomd.cluster_analysis.ClusterStructure()

    lin_protocol = espressomd.lees_edwards.LinearShear(
        initial_pos_offset=-0.02, time_0=0., shear_velocity=0.)
    system.lees_edwards.set_boundary_conditions(
        shear_direction="y", shear_plane_normal="x", protocol=lin_protocol)

    np.random.seed(1)

    def setup_polymer_through_shear_boundary(self, offset_factor):
        # place particles on a line across a Lees-Edwards periodic boundary
        for x in np.arange(-0.2, 0.21, 0.01)[::-1]:  # <--- different results without `[::-1]`
            x += 1e-7  # small offset due to Lees-Edwards edge case at x=0
            pos = [x, 1.1 * x, 1.2 * x]
            vel = [0., 0., 0.]
            # for particles beyond the shear boundary, place them initially in
            # the central box with an x-velocity that will cause them to cross
            # the shear boundary in 1 time step; the shear y-offset is exactly
            # compensated to form a straight line in unfolded coordinates
            # (only when `offset_factor = 1`)
            if x < 0.:
                pos[0] += 0.2
                vel[0] -= 0.2 / self.system.time_step
                pos[1] -= offset_factor * self.lin_protocol.initial_pos_offset
            self.system.part.add(pos=pos, v=vel)
        self.system.integrator.run(1)

    def test_lees_edwards(self):
        # Cluster analysis and fractal dimension require special treatment with
        # Lees-Edwards. The center of mass depends on whether the reference
        # particle (arbitrarily chosen) is positioned on the left or the right
        # of the shear plane, and the position offset is counted twice.

        # Place particles on a line (crossing Lees-Edwards periodic boundary)
        self.setup_polymer_through_shear_boundary(2.)  # <--- magic value

        dc = espressomd.pair_criteria.DistanceCriterion(cut_off=0.13)
        self.cs.pair_criterion = dc
        self.cs.run_for_all_pairs()
        self.assertEqual(len(self.cs.clusters), 1)

        c = list(self.cs.clusters)[0]
        # Discard cluster id
        c = c[1]

        # Center of mass should be at origin
        np.testing.assert_allclose(
            c.center_of_mass(), [1e-7, 1.1e-7, 1.2e-7], rtol=1E-8)

        # Longest distance
        ref_dist = self.system.distance(
            self.system.part.by_id(0),
            self.system.part.by_id(len(self.system.part) - 1))
        self.assertAlmostEqual(c.longest_distance(), ref_dist, delta=1E-8)

        # Radius of gyration
        rg = 0.
        com_particle = self.system.part.by_id(len(self.system.part) // 2)
        for p in c.particles():
            rg += self.system.distance(p, com_particle)**2
        rg /= len(self.system.part)
        rg = np.sqrt(rg)
        self.assertAlmostEqual(c.radius_of_gyration(), rg, delta=1E-12)

        self.system.part.clear()
        self.cs.clear()
        self.system.time = 0.

        # The fractal dimension of a line should be 1
        if espressomd.has_features("GSL"):
            self.setup_polymer_through_shear_boundary(np.log(3) / np.log(2))  # <--- magic value
            self.assertAlmostEqual(c.fractal_dimension(dr=0.)[0], 1, delta=0.05)


if __name__ == "__main__":
    ut.main()