Add short range neighbor methods

[jhossbach]

Fixes #4399, fixes #4438

Description of changes:

• CellStructure.cpp: Added the CellStructure::run_on_particle_short_range_neighbors() method which runs a kernel function over all particles inside the cell and its neighbors of a given particle
• cells.hpp: Added the mpi_get_short_range_neighbors() function to execute a parallel search
  • can be called from the Python interface via system.cell_system.get_neighbors.get_neighbors(p, distance)
• energy.hpp: Added the compute_non_bonded_pair_energy() function which returns both the short-range Coulomb interactions plus the non-bonded energy contributions of two particles
  • can be called from the Python interface via system.analysis.particle_energy(p)
• Reaction methods can use the new neighbor search algorithm with constructor argument search_algorithm="parallel" (default is search_algorithm="order_n"); on 1 MPI rank the original order N algorithm is faster since the parallel algorithm introduces some overhead due to the ghost update (this overhead is negligible with 2 or more MPI ranks)

[jhossbach]

The most recent version can now be found here: https://github.com/jhossbach/espresso/tree/short_range_neighbors_patch
@stekajack: Most tests in the unit test produce the expected results now, but there is still a problem with PBC neighbors which I will look at during the next coding day.

[jhossbach]

Left ToDo's :

• Write unit_tests/short_range_neighbors_test.cpp to test that neighbor cells are chosen correctly
• Add python interface for the short range neighbors (see Efficient search for interacting particles #4399)
• Add python interface for the short range energy contribution
• Test short range neighbors via python interface (or as a cpp unit test within in unit_tests/short_range_neighbors_test.cpp)
• Test particle_short_range_energy_contribution within python interface

[jhossbach]

One of the changes introduced an error with the mpi_... functions, I will check and push an updated version.

[jngrad]

Looks good to me! Once the remaining proposed changes are in, we can move forward with the merge.

[jhossbach]

There is a problem with getting the right distances when using more than one MPI rank, this needs to be adressed first

[jngrad]

@RudolfWeeber could you please have a look at this PR? Somehow running a short-range loop (e.g. to get the list of neighbors) updates the ghost information, but only in the cellsystem, so we end up with a corrupted particle node cache. This is why running analysis functions that act on particle positions after particles have been moved return results based on the old position, unless we run the integrator with 0 steps.

To reproduce the issue, run the test on 3 MPI ranks. This triggers a segfault due to a corrupted particle node cache.

When replacing

  if (pnode == this_node) {
    assert(cell_structure.get_local_particle(part));
    return *cell_structure.get_local_particle(part);
  }

by

  if (pnode == this_node) {
    auto const p_ptr = cell_structure.get_local_particle(part);
    if (p_ptr) {
      return *p_ptr;
    }
  }

in src/core/particle_data.cpp:get_particle_data(), the segfault moves to the particle insertion code. When commenting out the block

  if (cell_structure.get_resort_particles()) {
    cells_update_ghosts(global_ghost_flags());
  }

that was introduced by this PR in src/core/cells.cpp:mpi_get_short_range_neighbors_local(), the segfault persists.

The more fundamental issue here is that doing

system.part.all().pos = np.random.random((20, 3)) * system.box_l

does not update the particle properties, but instead queues the change until the next integration loop. So at the minimum one has to do

system.integrator.run(0, recalc_forces=False)

before attempting to update particle positions again. This also explains why the ParticleHandle property node still returns the old MPI rank after moving the particle to a different MPI domain, until the next integration.

[RudolfWeeber]

Please try on_observable_calc() from event.cpp rather than the manual resort call at the beginning of the neigbor search.

[RudolfWeeber]

The resort which was already present means that particles can move from one node to an other. Then the particle node cache needs to be invalidated, which is done in the above mentioned funciton.

[jhossbach]

The assertion in the check_neighbors function fails when running with 5 and 7 MPI cores. Also, we should implement a runtime check to ensure that the distance is less than the cell size, otherwise particles that are within distance but not in the neighborlist of that cell will not be included

[jhossbach]

I added a sanity check to prevent looking for particles for distances further than the cell size.
I also adapted the unit test, can you have another look @jngrad?

[jngrad]

The HybridDecomposition cell system is now partially supported thanks to @pkreissl. This was achieved by iterating over the N-square cells and injecting the local regular decomposition cell as a neighbor. This means that in a system with ions on a regular decomposition cell system and colloids on a N-square cell system, ions can detect any colloid without any distance restriction, however colloids can only detect ions in the regular cell of the same MPI node. I tried adding the regular cell neighbors on the local node to the local N-square cell, but then I get double counting since the regular cell neighbors already contain the N-square cells. To fix that, one would have to create temporary vectors in the hybrid decomposition constructor. However, there is another issue: since the colloid cell has no geometric boundary, a colloid fixed in space can potentially detect all ions in regular cell of MPI rank 0 at time step 0, and all ions in regular cell of MPI rank 1 at time step 1, if the N-square cells are resorted. To fully support HybridDecomposition, one probably would have to add all regular cells as neighbors of the N-square cells.