Cf JuliaMolSim activities - please join the slack if you haven’t yet!
1 Like
Since this is not an urgent transition fir me I will wait until the community has agreed on an abstract interface for Atoms and Structures at which point I will retire about half to 2/3 of JuLIP including the neighbour lists and much of the structure management and focus just in potentials again.
1 Like
@lmiq thanks for taking time to write this package.
Looks very interesting to test this using SPH too 
I have one question though, is the goal of the package you are making to be able to perform molecular simulations or is it simply to be able to generate neighbour lists?
Kind regards
1 Like
Actually the goal is the computation of the objective function and gradient for Packmol. That is closer to what is needed for a MD simulation than for neighbor lists. My comparison there is just because that was a good package to compare with.
Ps. What is SPH?
1 Like
Ah I see.
SPH stands for Smoothed Particle Hydrodynamics, which is a numerical particle simulation method. Can be used in different fields, I mostly use it in fluid dynamics. Has some similarities to molecular dynamics I suppose, but in SPH we solve the Navier-Stokes equation for each particle to model flow fields, pressure, forces on objects etc.
Nice visualization here
Kind regards
1 Like
I guess it could be used for that, if those simulations require the computation of interactions within a cutoff.
@lmiq
There is a similar package called CellLists.jl
https://github.com/jaantollander/CellLists.jl
It would be nice to know how CellListMap.jl is different from CellLists.jl in terms of performance !
The easiest comparison is probably with the neighbour lists functions of both packages:
julia> using CellListMap
julia> import CellLists
julia> n, d, r = 10_000, 3, 10.; p = 100*rand(n,d);
julia> using BenchmarkTools
julia> @btime let p = $p, r = $r
c = CellLists.CellList(p, r)
CellLists.near_neighbors(c,p,r)
end
114.365 ms (3387927 allocations: 361.06 MiB)
186283-element Vector{Tuple{Int64, Int64}}:
(2190, 3458)
(2190, 3516)
⋮
(9932, 2019)
(9932, 9024)
julia> x = [ p[i,:] for i in 1:n ];
julia> @btime CellListMap.neighbourlist($x,10,parallel=false)
11.073 ms (4979 allocations: 8.97 MiB)
186283-element Vector{Tuple{Int64, Int64, Float64}}:
(1, 383, 8.399067107896867)
(1, 1316, 9.936877451715944)
⋮
(9755, 7861, 8.58091039978657)
(9755, 8268, 9.898527922742463)
julia> @btime CellListMap.neighbourlist($x,10,parallel=true)
5.533 ms (5288 allocations: 19.61 MiB)
186283-element Vector{Tuple{Int64, Int64, Float64}}:
(1, 383, 8.399067107896867)
(1, 1316, 9.936877451715944)
⋮
(9755, 7861, 8.58091039978657)
(9755, 8268, 9.898527922742463)
julia> Threads.nthreads()
4
3 Likes
That sounds like a huge speed up with CellListMap.jl
Does CellListMap use bucketing approach also?? Means only searching for cells with indices with indices more than current index?
http://efekarakus.github.io/fixed-radius-near-neighbor/
I think this might give additional performance
2 Likes
Yes, it does for orthorhombic cells. For general triclinic cells it does not (yet?), it is a little bit trickier in that case. Indeed it can result in an almost 2x speedup (the inputs below are artificial in terms of the interface, just to force the cell type in the example):
julia> @btime let x = $x
box = Box(limits(x).limits .+ 10,10, UnitCellType=TriclinicCell)
cl = CellList(x,box,parallel=false)
CellListMap.neighbourlist(box,cl,parallel=false)
end;
17.324 ms (4976 allocations: 8.97 MiB)
julia> @btime let x = $x
box = Box(limits(x).limits .+ 10,10, UnitCellType=OrthorhombicCell)
cl = CellList(x,box,parallel=false)
CellListMap.neighbourlist(box,cl,parallel=false)
end;
9.880 ms (4976 allocations: 8.97 MiB)
(edit: that is not the only difference between the algorithms, but probably the most significant one. There are other small optimizations that can be done for orthorhombic cells).
1 Like
This is an update of the benchmark shown in the first post (which apparently cannot be edited anymore), for one recent version and also fixing the comparison with NearestNeighbours.jl, because the time of both packages if very much dependent on which set is chosen as the one for which the lists or trees will be built. In general building the list for the smallest set is the most efficient choice (which is made by default in the CellListMap.neighbourlist function).
Computing pairwise velocities of “galaxies”:
In-range neighbour search (the fastest of the two options of which set is used to build the list is shown. Parallelization with 4 physical cores - 8 threads; nl: NearestNeighbour.jl inrange function; cl: CellListMap.neighbourlist function):
----------------------------------
N1 = 1 ; N2 = 1000000 ; PASS TEST = true
nl (x,y): 102.679 ms (1000761 allocations: 83.96 MiB)
cl serial (y,x): 133.007 ms (39 allocations: 23.01 MiB)
cl parallel (y,x): 54.008 ms (98 allocations: 23.05 MiB)
----------------------------------
N1 = 10 ; N2 = 1000000 ; PASS TEST = true
nl (x,y): 146.351 ms (1006920 allocations: 84.24 MiB)
cl serial (y,x): 131.896 ms (68 allocations: 23.34 MiB)
cl parallel (y,x): 67.706 ms (136 allocations: 23.63 MiB)
----------------------------------
N1 = 100 ; N2 = 1000000 ; PASS TEST = true
nl (x,y): 294.152 ms (1064619 allocations: 86.89 MiB)
cl serial (y,x): 139.902 ms (288 allocations: 27.51 MiB)
cl parallel (y,x): 58.175 ms (369 allocations: 32.01 MiB)
----------------------------------
N1 = 1000 ; N2 = 1000000 ; PASS TEST = true
nl (x,y): 831.164 ms (1491040 allocations: 106.46 MiB)
cl serial (y,x): 259.863 ms (2216 allocations: 49.02 MiB)
cl parallel (y,x): 100.269 ms (2322 allocations: 71.60 MiB)
----------------------------------
N1 = 10000 ; N2 = 1000000 ; PASS TEST = true
nl (x,y): 2.318 s (3027684 allocations: 224.04 MiB)
cl serial (x,y): 761.103 ms (59024 allocations: 305.87 MiB)
cl parallel (x,y): 539.894 ms (59172 allocations: 454.62 MiB)
----------------------------------
N1 = 100000 ; N2 = 1000000 ; PASS TEST = true
nl (y,x): 6.452 s (987257 allocations: 1.51 GiB)
cl serial (x,y): 7.514 s (69905 allocations: 1.55 GiB)
cl parallel (x,y): 4.441 s (70059 allocations: 3.16 GiB)
----------------------------------
3 Likes
Out of curiosity, do you believe it would be hard to apply the same algorithm you are using now on GPU?
Kind regards
I think it is doable, and it is on my plans. But I don’t have much experience porting things to GPUs, so I cannot wright now anticipate how much work is needed to do that effectively. (MD codes run on GPUs well, so one way or the other that has to be possible)
I see thanks for the comment, very interesting. I once tried to write some simulation code with the idea of being able to use the same code for cpu and gpu, but I really struggled with it. Just wanted to hear about others experience.
Still, very impressive what CellListsMap is able to do right now!
Kind regards
1 Like
Here’s a snapshot of a simulation using Molly.jl with CellListMap.jl used to obtain neighbour lists.
It’s the designed protein Foldit1 (PDB ID 6MRR) in the a99SB-ILDN force field with explicit solvent (not shown).
11 Likes
Thanks for sharing! If you don’t mind I will add that to the showcase in CellListMap
One curiosity: do you have already a good PME implementation to deal with the long-range electrostatics? I searched some time ago but there didn’t seem to be a package available for that yet, and I was wandering if it was the case of implementing one.
Yes sounds good!
No not yet, I am using the reaction field approximation for this example. I think it’s a case of implementing it, see some discussion at Electrostatics in Julia.
1 Like
I am happy to announce that the article describing this package is now published, and available for free access for 50 days here:
CellListMap.jl: Efficient and customizable cell list implementation for calculation of pairwise particle properties within a cutof. Computer Physics Communications, Volume 279, October 2022, 108452.
The pre-print is also available at arXiv, here.
Associated to the paper, there are some Pluto notebooks explaining each of the examples shown at: GitHub - m3g/CellListMapArticleCodes: Code blocks of the CellListMap paper..
One more time, I feel obliged to acknowledge all the contributions from this forum (and from anonymous and not so anonymous referees of the article ), which allowed me to develop this and other Julia tools. It has been a very productive and exciting programming journey.
13 Likes
Congratulations! A very interesting and informative article to read.
1 Like