LB+particles: guard against missed coupling due to round-off error

[RudolfWeeber]

Fixes #4825

Description of changes:

• Bugfix: particles outside the simulation box are now properly coupled using PBC. The coordinates are folded before considering shifted positions in the LB particle coupling code. Also, a test that fails without the fix is added. To my understanding, if the particle position is folded, and then all combination of folded_pos[i] +/- box_l[i] for all Cartesian directions are considered, both the particle in the primary box and all potential halo regions are caught.

[RudolfWeeber]

in algorithsm/periodic_fold.hpp, first the lower, then the upper limit is checked, so here
-epsilon + box_l = box_l (due to round-off error)
but box_l >= box_l so
folded_pos = box_l - box_l =0

Let's consider the other case:
box_l +epislon is folded to
box_l +epsilon -box_l
even if this evaluates to 0 (rather than epsilon) due to rounding, the resulting coordinate is still folded, ie., 0<=x<box_l

box_l-epsilon is < box_l so the folding code does nothing, the coordinate is already folded.

So as I understand it, Algorithm::periodic_fold() is correct in all cases. Can someone please double check this?

[christophlohrmann]

I agree with the code changes and all arguments made.
Still, a cpp unit test for algorithm/periodic_fold that explicitly tests the cases mentioned in @RudolfWeeber's comment would support the argument better than my agreement.
Also the python test could look not only at the fact that the particle is coupled, but also could check that the particle is coupled correctly, i.e. momentum conserving.

[jngrad]

@RudolfWeeber's argument on the periodic folding w.r.t. roundoff errors is supported by the C++ unit tests:

espresso/src/core/unit_tests/periodic_fold_test.cpp

Lines 101 to 110 in 2d65667

	/* Corner left */
	{
	auto const x = std::nextafter(0., -1.);
	auto const box = 10.;
	auto const res = periodic_fold(x, 0, box);
	BOOST_CHECK(res.first >= 0.);
	BOOST_CHECK(res.first <= box);
	BOOST_CHECK(std::abs(res.first - x + res.second * box) <=
	4. * std::numeric_limits<double>::epsilon());
	}

Here res.first (position) is 0.0 and res.second (image box count) is 0, so the folding was applied twice. This test also checks for the case box+epsilon, and for a non-zero image box counter.