NeighbourLists

[cortner]

discuss cell list imlementation here vs NeighbourLists.jl - worth combinig code an efforts?

[m3g]

I have split my implementation into this package: http://github.com/m3g/CellListMap.jl

It only supports orthorhombic cells. Could we compare the performance of both to see if we are more or less on the same page? I have put some simple examples there, but without the documentation I am not sure how to compare to your implementation.

[cortner]

Sorry for the delay. CellListMap.jl is vastly more efficient ....

[ Info: Timings - CellListMap.jl:
        CellList:   6.109 ms (49 allocations: 4.24 MiB)
   map_pairwise!:   20.064 ms (321 allocations: 11.47 KiB)
[ Info: Timings - NeighbourLists.jl:
  NeighbourLists:   1.091 s (800095 allocations: 293.98 MiB)
  pairs iterator:   11.660 ms (0 allocations: 0 bytes)

, so much so that I wonder whether I messed something up at some point since when I first wrote it, it was about on par with a very popular (and supposedly fast) neighbourlist written in C++. But who knows; I'll do a bit more debugging to see if I can find an obvious performance bottleneck.

Couple options

• if your CellListMap.jl can eventually support arbitrary cell shapes and PBC on or off, then I'd be thrilled to retire NeighbourLists.jl.
• The alternative might be to figure out what you did right and I did wrong and adapt NeighbourLists.jl accordingly

[cortner]

after a bit of profiling it appears that ALL the work is done when converting the cell-list to a pair list. Since you do that implicitly in the map_pairwise! function you are not really avoiding it either... So I wonder whether the way you choose your cells and cell sizes is just vastly superior. If I remember correctly, my cell size is the cutoff radius. And it is quite possible that this is an ok heuristic for molecular simulation but that it is terrible for random datapoints.

[cortner]

is there a way to move this issue over to CellListMap.jl ?

[lmiq]

NeighbourLists: 1.091 s (800095 allocations: 293.98 MiB)

this kind of allocations suggest that there is some type instability in the code, which may be causing most of the slowdown.

Apart from that, I have tweaked quite a few things. I run only over the cells which are "forward" on the list to avoid repetitions (thus over 13 cells per cell, not 26). The cell size is an adjustable parameter , but by default this parameter is 1 and the cell size is the same as the cutoff. There are many other small optimizations that I have done in the last few weeks also.

I will effectively implement support for general periodic boundary conditions, I expect very soon, this is on my priority list. Also, I will try to implement the ideas on this paper as well: https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiClvCWqLrxAhWMpJUCHZZ8D-wQFjABegQIAxAD&url=https%3A%2F%2Farxiv.org%2Fabs%2F1804.06231&usg=AOvVaw2m5jTOkSYfeTx5pLaehCar

But my goal is not to provide the best strategy for MD simulations. Actually for MD simulations the best approach, the Verlet Lists, is easier: build a list of neighbours in the first step (maybe using cell lists), and update these lists after some steps only for the particles that moved outside the boundaries. This is much faster than the cell lists approach, if you can guarantee that there is correlation between the positions of the particles among successive frames. This is not my case though, as for analysis routines the frames are to far from each other in time, and for Packmol, there might be no correlation between the positions of the particles in different steps of the optimization procedure.

[cortner]

We never bothered with verlet lists since we are interested in potentials that are vastly more expensive than a neighbourlist assembly.

that’s also why I never noticed that there might be a problem :) - or indeed if I had might not have bothered to fix it :). But now I’m curious

[cortner]

The profiling doesn’t show an allocation where the bottleneck is. I thought this meant there is no type instability. Otherwise I 100% agree it looks just like that

[lmiq]

That is strange. Maybe a good pass of @code_warntype in different functions finds something.

[cortner]

I will definitely do that next once I find the time :). Atm we have school closure due to heatwave so no chance…

[lmiq]

@cortner, just to inform that version 0.5.1 now supports arbitrary periodic boundary conditions (and should be somewhat faster yet, particularly for dense systems).

[cortner]

that's really great, thanks for mentioning it. I'll run it through my test set and maybe I can retire my NeighbourLists.jl package.

[lmiq]

I am working actively on the package, and it will become faster still (it will fall back to different strategies of computation depending on system size and density). I still have to set a scheme for if someone does not want periodic boundary conditions at all, but that will be easy.

[m3g]

Not sure if you had time to test anything, but I noticed that the general PBCs where failing in many cases. Now I uploaded a new version, 0.5.2, which will raise errors if the cell cannot be handled. In particular I can only handle correctly cells with all positive lattice vectors and which are not collapsed relative to the cutoff. At some point I will try to follow the same rules either Gromacs or Lammps follow, but that is much more cumbersome than I initially imagined.

[cortner]

Sorry I didn’t get to it yet, too many other things going on. But I will for sure look into it.

[lmiq]

The package provides a neighbourlist function that implements this calculation. Without periodic boundary conditions, just do:

julia> x = [ rand(3) for _ in 1:10_000 ];

julia> CellListMap.neighbourlist(x,0.05)
24778-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 62, 0.028481068525796384)
 ⋮
 (9954, 1749, 0.04887502372299809)
 (9974, 124, 0.040110356034451795)

or CellListMap.neighbourlist(x,y,r) for computing the lists of pairs of two sets closer than r.

The returning array contains tuples with the index of the particle in the first vector, the index of the particle in the second vector, and their distance.

If periodic boundary conditions are used, the Box and CellList must be constructed in advance:

julia> x = [ rand(3) for _ in 1:10_000 ]; 

julia> box = Box([1,1,1],0.1);

julia> cl = CellList(x,box);

julia> CellListMap.neighbourlist(box,cl)

julia> CellListMap.neighbourlist(box,cl)
209506-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 121, 0.05553035041478053)
 (1, 1589, 0.051415489701932444)
 ⋮
 (7469, 7946, 0.09760096646331885)

(this is part of the documentation now, and will be released exactly like this in v0.7)