m3g / celllistmap.jl Goto Github PK

See: JuliaDynamics/Agents.jl#659

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

We are currently looking at implementing replica exchange methods in Molly (JuliaMolSim/Molly.jl#64), which involves running a number of simulations in parallel. In the case of 4 replicas you might want to use 4 cores per replica if you had 16 cores, with the neighbour lists for each replica being separate.

It seems in CellListMap.jl that you can set parallel=true and this will use a number of cores up to Threads.nthreads():

Is there a way to say not just parallel=true but to set a maximum number of threads? If not, how easily could this feature be added? This would make it easier to use CellListMap.jl with multiple replica simulations.

Hi Again, one more question:

particles in biology are organisms with varying traits. This means for example, that in a spatial food web simulation the distance at which an organisms reacts to another organism is dependent on .e.g. the organism size and can vary very much.

If I run a model with CellListMaps I now need to put the maximum reaction distance of all organisms as cutoff and then implement a per organism condition in the function that gets called by map_pairwise! to account for different reaction distances.

I see that cutoff in e.g. PeriodicSystem is of type number. But would it be possible to optionally provide an array of per particle reaction distances that sets the cutoff per organism/particle?

UpdateCellList! needs more extensive testing with box size changes and coordinate changes, in parallel or not.

This can be useful for a self-avoiding random walk simulation, for example.

See here: https://discourse.julialang.org/t/celllistmap-jl-is-it-suitable-for-detecting-collision-of-particles/

Assume I simulate individual organisms. Sometimes they die and their mass goes to zero. In that case I don't want the effect of them to be calculated. I now do that by

function f(i,j,d2,d)
	d[i]+= u[i]>0.0 ? exp(-d2*dc)*u[j]*(u[j]/u[i])^ac : 0.0
	return d
end

in the function called by CellListMap.map_pairwise! which speeds it up somewhat. But can I also keep a list of organisms which are dead and make CellListMap.map_pairwise! just skip the call to the function altogether? When I have a birth of an organism, I take one dead organism and revive it with the new genetic makeup (traits) by setting its mass>0.

So the short question is, can I somehow keep a list of which organisms are alive without having to redo the entire SArray of positions for every death/birth event?

Hi, I have been using CellListMap for a while, and a few days ago updated the package. To my surprise, now disperse systems as simple as 1 or 2 coordinates throw an error that didn't happen before. Here is son simple code that can reproduce the error:

r = [[1.0,1.0,1.0]]
system = InPlaceNeighborList(x=r, cutoff=3.0, parallel=false)
list = neighborlist!(system)
update!(system, r)

This also happens when data is a bit dispersed. For instance, with these coordinates it happens too: r = [[1.0,1.0,1.0], [10.0, 1.0, 1.0], [3.0, 1.0, 1.0]]

The error is the following:

UNIT CELL CHECK FAILED: distance between cell planes too small relative to cutoff.
ArgumentError:  Unit cell matrix does not satisfy required conditions.

It is weird because list = neighborlist!(system) is working fine.

For now, I will revert to the old version, but this is something that should be fixed.

Best regards.

Hi, I wanted to try getting some gradients from a function involving map_pairwise as I saw on the docs that automatic differentiation was available. The issue is my input consists in hundreds of thousands of variables (a 3x~1e5 Matrix) and my output is a loss score, so ForwardDiff is quite inefficient. I tried just replacing it with ReverseDiff but I got this

ERROR: LoadError: ArgumentError: cannot reinterpret `ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}` as
 `SVector{3, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}`, type `SVector{3, ReverseDiff.TrackedReal
{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}` is not a bits type 

and with Zygote

ERROR: LoadError: Compiling Tuple{CellListMap.var"##set_number_of_batches!#38", Bool, typeof(CellListMap.set_number_of_batches!), CellList{3, Float64}, Tuple{Int64, Int64}}: try/c
atch is not supported.

My question would be first, if it is even possible to use reverse-mode differentiation with CellListMap? If so, is it possible to add some examples to the docs on how to do so? The type-conversion trick used for ForwardDiff does not work.
Thanks in advance for your help.

With Julia v1.7.2 and CellListMap v0.7.20, the example from the docs works fine:

using CellListMap, Unitful, StaticArrays
cutoff = 0.1u"nm"
box = Box([1.0, 1.0, 1.0]u"nm", cutoff)
x = [rand(typeof(cutoff), 3) for _ in 1:1000]
cl = CellList(x, box; parallel=false)

However when trying it with a triclinic box:

box2 = Box([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]u"nm", cutoff)
cl2 = CellList(x, box2; parallel=false)

ERROR: LoadError: DimensionError: nm^-2 and 0.8425463402211867 are not dimensionally compatible.
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐋^-2, Unitful.FreeUnits{(nm^-2,), 𝐋^-2, nothing}}}, x::Float64)
    @ Unitful ~/.julia/packages/Unitful/SUQzL/src/conversion.jl:112
  [2] macro expansion
    @ ~/.julia/packages/StaticArraysCore/gkLqH/src/StaticArraysCore.jl:79 [inlined]
  [3] convert_ntuple
    @ ~/.julia/packages/StaticArraysCore/gkLqH/src/StaticArraysCore.jl:75 [inlined]
  [4] SVector{3, Quantity{Float64, 𝐋^-2, Unitful.FreeUnits{(nm^-2,), 𝐋^-2, nothing}}}(x::Tuple{Float64, Float64, Float64})
    @ StaticArraysCore ~/.julia/packages/StaticArraysCore/gkLqH/src/StaticArraysCore.jl:111
  [5] StaticArray
    @ ~/.julia/packages/StaticArrays/8Dz3j/src/convert.jl:160 [inlined]
  [6] _solve
    @ ~/.julia/packages/StaticArrays/8Dz3j/src/solve.jl:17 [inlined]
  [7] \
    @ ~/.julia/packages/StaticArrays/8Dz3j/src/solve.jl:1 [inlined]
  [8] wrap_cell_fraction
    @ ~/.julia/packages/CellListMap/zdcxg/src/CellOperations.jl:46 [inlined]
  [9] wrap_to_first
    @ ~/.julia/packages/CellListMap/zdcxg/src/CellOperations.jl:80 [inlined]
 [10] wrap_to_first
    @ ~/.julia/packages/CellListMap/zdcxg/src/CellOperations.jl:98 [inlined]
 [11] add_particles!(x::Vector{Vector{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}, box::Box{TriclinicCell, 3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(nm^2,), 𝐋^2, nothing}}, 9}, ishift::Int64, cl::CellList{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}})
    @ CellListMap ~/.julia/packages/CellListMap/zdcxg/src/CellLists.jl:871
 [12] UpdateCellList!(x::Vector{Vector{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}, box::Box{TriclinicCell, 3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(nm^2,), 𝐋^2, nothing}}, 9}, cl::CellList{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}, aux::Nothing; parallel::Bool)
    @ CellListMap ~/.julia/packages/CellListMap/zdcxg/src/CellLists.jl:787
 [13] #UpdateCellList!#49
    @ ~/.julia/packages/CellListMap/zdcxg/src/CellLists.jl:679 [inlined]
 [14] CellList(x::Vector{Vector{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}, box::Box{TriclinicCell, 3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(nm^2,), 𝐋^2, nothing}}, 9}; parallel::Bool, nbatches::Tuple{Int64, Int64})
    @ CellListMap ~/.julia/packages/CellListMap/zdcxg/src/CellLists.jl:495
 [15] top-level scope
    @ ~/dms/molly_dev/cnf_err.jl:11
in expression starting at /home/jgreener/dms/molly_dev/cnf_err.jl:11

Stripping the units in the triclinic case also works fine.

Currently, UpdateCellList! clears all cells and adds particles again. For use cases like SPH simulations where only small changes in position happen between updates, it will be much faster to only update cells that changed, using the existing cell lists from the previous step.

Remember to add testing functions for the case the inputs are given as matrices. And add the possibility of this input to the documentation.

Hello, thanks for the great package!

This is my second attempt at building something with this package.

I have a very simple code for doing Brownian Dynamics using a continuous pseudo hard-sphere potential in 3D.
This is the code in a gist

However, when the energy and virial get updated in the energy_and_forces! function, which is summarized here

function energy_and_forces!(x,y,i,j,d2,cutoff,output::EnergyAndForces)
    d = sqrt(d2)
    r = y - x
    if d < cutoff
        (uij, fij) = pseudohs(d)
        sums = @. fij * r / d
        output.virial += dot(sums, r)
    end
    output.energy += uij
    output.forces[i] = @. output.forces[i] + (fij * r / d)
    output.forces[j] = @. output.forces[j] - (fij * r / d)

    return output
end

I always get very large values for energies and the virial. I have been trying to understand what the problem is here, but I don't know where to look into further. Any help is extremely appreciated.

Hi, I am using CellListMap to apply a couple of functions to a catalog of ~16M objects. The general idea is as follows:

using Random
using CellListMap
mutable struct Catalog{T}
    pos::AbstractMatrix{T}
    vel::AbstractMatrix{T}
    r_min::AbstractVector{T}
    param_1
    param_2
end #struct


function update!(sat_pos, cen_pos, sat_id, cen_id, d2, catalog::Catalog{T}) where T
    param_1 = catalog.param_1
    param_2 = catalog.param_2

    catalog.pos[1, sat_id] = cen_pos[1] + randn(param_1 * param_2)
    catalog.pos[2, sat_id] = cen_pos[2] + randn(param_1 * param_2)
    catalog.pos[3, sat_id] = cen_pos[3] + randn(param_1 * param_2)

    for ax = 1:3
        catalog.vel[ax, sat_id] = catalog.vel[ax, sat_id] + param_2
    end #for
    catalog
end #func


function inner_1!(cen_pos, sat_pos, cen_id, sat_id, d2, catalog::Catalog{T}) where T<:Real
            catalog = update!(sat_pos, cen_pos, sat_id, cen_id, d2, catalog::Catalog{T})
            catalog.r_min[sat_id] = d2
    catalog
end #func

function inner_2!(cen_pos, sat_pos, cen_id, sat_id, d2, catalog::Catalog{T}) where T<:Real
            catalog = update!(sat_pos, cen_pos, sat_id, cen_id, d2, catalog)
            catalog.r_min[sat_id] = d2
    catalog
end #func


function apply(catalog::Catalog{T},
        params::AbstractVector,
        box_size::AbstractVector;
        show_progress = true) where T
    #catalog = copy_output(catalog)
    box = Box(box_size, params[1])
    cl = CellList(catalog.pos, catalog.pos, box)
    catalog.param_1 = params[3]
    catalog.param_2 = params[end]
    CellListMap.map_pairwise!((x, y, i, j, d2, catalog) -> inner_1!(x, y, i, j, d2, catalog),
                              catalog, box, cl, parallel=false, show_progress=show_progress)
   # UP TO THIS POINT THE CODE RUNS EXTREME ALLOCATION SEEMS TO BE FROM HERE ON
    box = Box(box_size, params[2])
    cl = CellList(catalog.pos, catalog.pos, box)
    catalog.param_1 = params[4]
    CellListMap.map_pairwise!((x, y, i, j, d2, catalog) -> inner_2!(x, y, i, j, d2, catalog),
                               catalog, box, cl, parallel=false, show_progress=show_progress)

    catalog

end #func

function main()
    Nobjects = 16000000
    cat = Catalog(2000f0 * rand(Float32, (3, Nobjects)), 2000f0 * randn(Float32, (3, Nobjects)), 2000f0 * rand(Float32, Nobjects), 0., 0.)
    box_size = [2000. for _ = 1:3]
    params = [1., 1., 1., 1., 1.]
    apply(cat::Catalog{Float32},
        params::AbstractVector,
        box_size::AbstractVector;
        show_progress = true)

end #func

main()

The catalog on disk occupies 1.4 Gb and while I use other large arrays, I wouldn't expect the code to need >300Gb. And it seems that it fails when updating the cell list, I just get Killed most of the time

Hi!

Thank you for this incredible tool for working with particle based cases. I am developing a post-processing tool https://github.com/AhmedSalih3d/PostSPH.jl directed at a software package called DualSPHysics (https://dual.sphysics.org/).

I have been using the current implementation of CellListMap (CPU based) and found that it is much faster than anything I have been able to write from hand, so I would of course want to make use of it in my package when post-processing different quantities.

The reason I ask about GPU is that it would be quite impressive being able to switch between CPU/GPU based on task and in the future I hope to write some small SPH simulations, which would greatly benefit from a GPU implementation of this library.

Please do not take me making this issue as a "demand", I just share in your passion about efficient neighbour algorithms for particle simulations! 😊

When systems are highly non-homogeneous, the number of cells can explode to be much greater than the number of particles:

using StaticArrays
using CellListMap

function points_in_sphere(r,N)
    p = SVector{3,Float64}[]
    for ip in 1:N
        θ = π*rand()
        ϕ = 2π*rand()
        x = r * sin(ϕ) * cos(θ)
        y = r * sin(ϕ) * sin(θ) 
        z = r * cos(ϕ)
        push!(p, SVector(x,y,z))
    end
    return p
end

julia> p = points_in_sphere(500, 100);

julia> box = Box(limits(p), 0.02)
Box{NonPeriodicCell, 3}
  unit cell matrix = [ 990.57, 0.0, 0.0; 0.0, 985.31, 0.0; 0.0, 0.0, 999.88 ]
  cutoff = 0.02
  number of computing cells on each dimension = [49530, 49267, 49995]
  computing cell sizes = [0.02, 0.02, 0.02] (lcell: 1)
  Total number of cells = 121997524527450

This makes it prohibitive to work with this kind of system. Can we find a way to not allocate empty cells at all? The full cell allocation would need to be lazy, since the construction of the box does not involve the actual coordinates of the system.

👋 Hi there!

As you may be aware, some of us over at JuliaMolSim have been developing AtomsBase, an interface package for specifying atomic geometries. This package seems like a phenomenal candidate to support it so that developers of packages making use of this tool can easily "hook it up" to their system!

I thought I'd put this on your radar, and would be happy to discuss the best way to go about implementing the interface, and/or contribute to a PR to do so!

In the chapter of periodic systems of the document it is said that update cell list could be disabled. I tried it on the non-periodic systems, and there is no 'update_lists' in the variables for map_pairwise!().
Attached is my code

using CellListMap
using StaticArrays
using BenchmarkTools
@inline function my_norm(v::AbstractVector{T}) where {T}
    n = zero(T)
    @simd for x in v
        n += x ^ 2
    end
    n = sqrt(n)
    return n
end
function cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage)
    d_new = my_norm(x + u_points[i] - y -u_points[j])
    λ = d_new - sqrt(d2)
    pair_dam = 0.0
    if λ > λᵦ
        pair_dam = 1.0 #- exp(-γ * (λ - λᵦ))
    end
    damage[i] += pair_dam * points_vol[j]
    damage[j] += pair_dam * points_vol[i]
    return damage
end
function map_damage!(damage_map,box,cl,u_points,points_vol,γ,λᵦ)
    damage_map .= 0.0
    map_pairwise!((x,y,i,j,d2,damage_map) -> cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage_map),damage_map,box,cl; update_lists=false)
end
N = 300000
ℓ = 0.05
γ = 1.0e4
λᵦ = 1.0e-6
points = rand(SVector{2,Float64},N)
box = Box(limits(points),ℓ)
cl = CellList(points,box)
u_points = 1e-4 * rand(SVector{2,Float64},N)
points_vol = rand(N)
damage_map = zeros(Float64,N)
map_pairwise!((x,y,i,j,d2,damage_map) -> cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage_map),damage_map,box,cl)
@btime map_damage!(damage_map,box,cl,u_points,points_vol,γ,λᵦ)

And there is the error information:

ERROR: MethodError: no method matching map_pairwise!(::var"#5#6"{Vector{SVector{2, Float64}}, Vector{Float64}, Float64, Float64}, ::Vector{Float64}, ::Box{NonPeriodicCell, 2, Float64, Float64, 4, Float64}, ::CellList{2, Float64}; update_lists=false)
Closest candidates are:
  map_pairwise!(::F, ::Any, ::Box, ::CellList; parallel, output_threaded, reduce, show_progress) where F at C:\Users\rydmaize\.julia\packages\CellListMap\BemzV\src\CellListMap.jl:90 got unsupported keyword argument "update_lists"
  map_pairwise!(::F1, ::Any, ::Box, ::CellListMap.CellListPair{V, N, T, Swap}; parallel, output_threaded, reduce, show_progress) where {F1, V, N, T, Swap} at C:\Users\rydmaize\.julia\packages\CellListMap\BemzV\src\CellListMap.jl:116 got unsupported keyword argument "update_lists"
Stacktrace:
  [1] kwerr(::NamedTuple{(:update_lists,), Tuple{Bool}}, ::Function, ::Function, ::Vector{Float64}, ::Box{NonPeriodicCell, 2, Float64, Float64, 4, Float64}, ::CellList{2, Float64})
    @ Base .\error.jl:165
  [2] map_damage!(damage_map::Vector{Float64}, box::Box{NonPeriodicCell, 2, Float64, Float64, 4, Float64}, cl::CellList{2, Float64}, u_points::Vector{SVector{2, Float64}}, points_vol::Vector{Float64}, γ::Float64, λᵦ::Float64)
    @ Main d:\Code\CellListMap\Update_lists.jl:26
  [3] var"##core#321"()
    @ Main C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:489
  [4] var"##sample#322"(::Tuple{}, __params::BenchmarkTools.Parameters)
    @ Main C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:495
  [5] _run(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; verbose::Bool, pad::String, kwargs::Base.Pairs{Symbol, Integer, NTuple{4, Symbol}, NamedTuple{(:samples, :evals, :gctrial, :gcsample), Tuple{Int64, Int64, Bool, Bool}}})
    @ BenchmarkTools C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:99
  [6] #invokelatest#2
    @ .\essentials.jl:731 [inlined]
  [7] #run_result#45
    @ C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:34 [inlined]
  [8] run(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; progressid::Nothing, nleaves::Float64, ndone::Float64, kwargs::Base.Pairs{Symbol, Integer, NTuple{5, Symbol}, NamedTuple{(:verbose, :samples, :evals, :gctrial, :gcsample), Tuple{Bool, Int64, Int64, Bool, Bool}}})
    @ BenchmarkTools C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:117
  [9] #warmup#54
    @ C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:169 [inlined]
 [10] warmup(item::BenchmarkTools.Benchmark)
    @ BenchmarkTools C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:168
 [11] top-level scope
    @ C:\Users\rydmaize\.julia\packages\BenchmarkTools\0owsb\src\execution.jl:575

Is it possible to make this option also possible for non-periodic systems? Thanks for your help!

Hi,

I've noticed that sometimes, the results of neighborlist are non-unique:

using StaticArrays
using CellListMap

l = SVector{3, Float32}[[-1.45, 0.0, 0.0], [0.0, 0.0, 0.0], [0.495, 0.0, 1.437], [1.235, -0.911, 1.838], [-1.788, 0.918, 0.25], [1.642, 1.124, -0.864], [-1.801, -0.26, -0.911], [0.154, 1.136, -1.827], [-1.749, -0.665, 0.704], [0.258, 2.118, -0.351], [0.554, 1.175, -0.813], [0.341, -0.928, -0.46], [0.133, 0.98, 2.268], [0.605, 0.98, 3.639], [0.414, -0.408, 4.228], [-0.476, 1.73, 1.938], [0.325, 2.106, 5.472], [0.065, 3.039, 3.988], [-1.16, 1.868, 4.519], [-0.089, 2.07, 4.463], [1.676, 1.185, 3.63], [0.531, -1.358, 3.421], [0.462, -0.512, 5.473]]

nl1 = neighborlist(l, 10.0) # returns a unique 253-element list
nl2 = neighborlist(l, 7.0) # returns a non-unique 260-element list; unique(nl2) has 250 elements

I want to calculate a damage variable which is somehow similar to the force in MD simulation with CellListMap.jl. In my problem, the pairwise damage is related to the deformation history of the pair, and I would like to store a historical value for each pair during calculation. Is there a possible way to do this with CellListMap.jl?
I will explain it with code. Attached is my code.

using CellListMap
using StaticArrays
using BenchmarkTools
@inline function my_norm(v::AbstractVector{T}) where {T}
    n = zero(T)
    @simd for x in v
        n += x ^ 2
    end
    n = sqrt(n)
    return n
end
function cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage)
    d_new = my_norm(x + u_points[i] - y -u_points[j])
    λ = d_new - sqrt(d2)
    pair_dam = 0.0
    if λ > λᵦ
        pair_dam = 1.0
    end
    damage[i] += pair_dam * points_vol[j]
    damage[j] += pair_dam * points_vol[i]
    return damage
end
function map_damage!(damage_map,box,cl,u_points,points_vol,γ,λᵦ)
    damage_map .= 0.0
    map_pairwise!((x,y,i,j,d2,damage_map) -> cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage_map),damage_map,box,cl)
end
N = 300000
ℓ = 0.05
γ = 1.0e4
λᵦ = 1.0e-6
points = rand(SVector{2,Float64},N)
box = Box(limits(points),ℓ)
cl = CellList(points,box)
u_points = 1e-4 * rand(SVector{2,Float64},N)
points_vol = rand(N)
damage_map = zeros(Float64,N)
map_pairwise!((x,y,i,j,d2,damage_map) -> cal_damage!(x,y,i,j,d2,u_points,points_vol,γ,λᵦ,damage_map),damage_map,box,cl)
@btime map_damage!(damage_map,box,cl,u_points,points_vol,γ,λᵦ)

In the function cal_damage!,

if λ > λᵦ
        pair_dam = 1.0
end

I would like to replace the λᵦ with a historical quantity, e.g., the maximum of λ in previous calculation, like

if λ > maximum of λ in [0,t]
        pair_dam = 1.0 
end

Is is possible with CellListMap.jl ? Thanks for your help.

See:

https://discourse.julialang.org/t/juliacall-pass-numpy-array-to-julia-function-as-vector-float64/87728/6?u=lmiq

Hi, I've been using CellListMap to compute some quick two-point functions (2pf). At first it seemed to work well but after some processing of my data, trying to recompute the 2pf hangs at about 70% eternally (apparently, at least more than one night). This is odd since before the processing, for the same number of points the 2pf calculation takes about 3s. After I Ctrl+C the calculation, it seems all Julia is stuck since everything is stuck until I relaunch the REPL.
Here I add the relevant parts of my code:

using CellListMap
using StaticArrays
using NPZ

function build_histogram!(d2,hist, bin_edges)
    ibin = searchsortedlast(bin_edges, sqrt(d2))
    #ibin = floor(Int,d) + 1
    if (ibin > 0) && ibin <= length(bin_edges)
        hist[ibin] += 1
    end #if
    return hist
end;

function compute_2pcf(catalog, box_size, bin_edges)
    sides = [eltype(catalog)(b) for b in box_size]
    #cutoff = box_size[1] / size(counts_ref, 1)
    cutoff = eltype(catalog)(bin_edges[end])
    box = Box(sides, cutoff)
    
    cl = CellList(catalog, box)
    
    # Initialize (and preallocate) the histogram
    hist = zeros(Int,size(bin_edges,1)-1);
    println("Counting pairs...")
    # Run calculation
    map_pairwise!(
        (x,y,i,j,d2,hist) -> build_histogram!(d2,hist, bin_edges),
        hist,box,cl; show_progress = true
    )
    println("Done")
    #hist = hist  / (length(satellites[1]) * (length(centrals[1])))
    N = size(catalog,1)
    hist = hist  / (N * (N - 1))
    norm = @. (4/3) * π * (bin_edges[2:end]^3 -bin_edges[1:end-1]^3) / (box_size[1] * box_size[2] * box_size[3])
    hist ./= norm
    centers = bin_edges[2:end] - bin_edges[1:end-1]
    centers, hist
end #

compute_2pcf(x, y, z, box_size, bin_edges) = compute_2pcf(Matrix(hcat(x,y,z)'), box_size, bin_edges)
    
box_size = @SVector [2000f0 for _ in 1:3]
bin_edges = 10 .^range(-2, stop = log10(50), length=11)
before  = npzread("before.npy")
after  = npzread("after.npy")

s, xi = compute_2pcf(before, box_size, bin_edges)                                                                                                                           
s, xi = compute_2pcf(after, box_size, bin_edges)   

output:

Counting pairs...                                                                                                                                                                  
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:06         
Done                                                                                                                                                                               
([0.013436729115920996, 0.03149129804938488, 0.07380530218921652, 0.17297548747273575, 0.4053979643592893, 0.9501202274834459, 2.226771039908675, 5.218822966551732, 12.23121401710
201, 28.665964967767582], [0.0, 1.6660056393630336e13, 2.038290235445229e13, 1.9276573437228812e13, 1.88991060147404e13, 1.8194857516864344e13, 1.4969815963181014e13, 7.5421693681
98844e12, 3.7787949312305874e12, 3.0175790787229033e12])                                                                                                                           
                                                                                                                                                                                   

Counting pairs...                                                                                                                                                                  
Progress:  71%|███████████████████████████████████████████████████████████████████████████████████████████████████▊                                        |  ETA: 0:00:00^CERROR: 
InterruptException:                                                                                                                                                                
Stacktrace:                                                                                                                                                                        
  [1] poptask(W::Base.IntrusiveLinkedListSynchronized{Task})                                                                                                                       
   @ Base ./task.jl:974                                                                                                                                                           
  [2] wait()
    @ Base ./task.jl:983
  [3] wait(c::Base.GenericCondition{Base.Threads.SpinLock}; first::Bool)
    @ Base ./condition.jl:130
  [4] wait(c::Base.GenericCondition{Base.Threads.SpinLock})
    @ Base ./condition.jl:125
  [5] _wait(t::Task)
    @ Base ./task.jl:308
  [6] sync_end(c::Channel{Any})
    @ Base ./task.jl:404 [7] macro expansion                                                                                                                                                     [35/1946]
    @ ./task.jl:477 [inlined]
  [8] map_pairwise_parallel!(f::var"#4#6"{Vector{Float64}}, output::Vector{Int64}, box::Box{OrthorhombicCell, 3, Float64, Float64, 9, Float64}, cl::CellList{3, Float64}; output_th
readed::Nothing, reduce::typeof(CellListMap.reduce), show_progress::Bool)
    @ CellListMap ~/.julia/packages/CellListMap/CQJ0D/src/CoreComputing.jl:148
  [9] map_pairwise_parallel!
    @ ~/.julia/packages/CellListMap/CQJ0D/src/CoreComputing.jl:135 [inlined]
 [10] #map_pairwise!#109
    @ ~/.julia/packages/CellListMap/CQJ0D/src/CellListMap.jl:98 [inlined]
 [11] compute_2pcf(catalog::Matrix{Float64}, box_size::SVector{3, Float32}, bin_edges::Vector{Float64})
    @ Main ./REPL[6]:13
 [12] top-level scope
    @ REPL[18]:1

I've checked there is no NaN or infinities in the after array.
I've also removed threads and it didn't even show the progress bar. When sending Ctrl+C nothing happened until I had to kill the whole process so no traceback there.

Thanks in advance for your help!

I have data from here AhmedSalih3d
repo .

And I found different results with brute force test.

path    = dirname(@__FILE__)
DF_FLUID = CSV.read(joinpath(path, "./FluidPoints_Dp0.02.csv"), DataFrame)
DF_BOUND = CSV.read(joinpath(path, "./BoundaryPoints_Dp0.02.csv"), DataFrame)
P1F = DF_FLUID[!,"Points:0"]
P2F = DF_FLUID[!,"Points:1"]
P3F = DF_FLUID[!,"Points:2"]
P1B = DF_BOUND[!,"Points:0"]
P2B = DF_BOUND[!,"Points:1"]
P3B = DF_BOUND[!,"Points:2"]
points = Vector{SVector{3,Float64}}()
for i = 1:length(P1F)
    push!(points,SVector(P1F[i],P3F[i],P2F[i]))
end
for i = 1:length(P1B)
    push!(points,SVector(P1B[i],P3B[i],P2B[i]))
end

H = 0.06788225099390856
system  = InPlaceNeighborList(x=points, cutoff = H, parallel=true)
update!(system, points)
list = neighborlist!(system)

101128 pairs

system  = InPlaceNeighborList(x=points, cutoff = H, parallel=false)
update!(system, points)
list = neighborlist!(system)

101128 pairs

 @inline function distance_condition(p1, p2, r) 
    sum(@. (p1 - p2)^2) ≤ r^2
 end

function brute_force(p, r)
    ps = Vector{Tuple{Int, Int}}()
    n = length(p)
    for i in 1:(n-1)
        for j in (i+1):n
            if @inbounds distance_condition(p[i], p[j], r)
                    push!(ps, (i, j))
            end
        end
    end
    return ps
end

brute_force(points, H)

98585 pairs

But if change a little bit:

system  = InPlaceNeighborList(x=points, cutoff = H - 0.0000001, parallel=false)
update!(system, points)
list = neighborlist!(system)

What can I do to get more reproducible results?

Also:

brute_force(points, H + sqrt(eps()))

98585 pairs

bf = brute_force(points, H + 0.00422)

98585 pairs

bf = brute_force(points, H + 0.00423)

118425 pairs

Something like this doesn't work, for no other reason than too specific dispatching:

julia> using CellListMap, StaticArrays

julia> x = [[0.0, 0.0]]; y = [SVector(0.0, 0.0)];

julia> limits(x, y)
ERROR: MethodError: no method matching limits(::Vector{Vector{Float64}}, ::Vector{SVector{2, Float64}})

Closest candidates are:
  limits(::AbstractVector{<:AbstractVector})
   @ CellListMap ~/.julia/dev/CellListMap/src/CellOperations.jl:296
  limits(::T, ::T) where T<:(AbstractVector{<:AbstractVector})
   @ CellListMap ~/.julia/dev/CellListMap/src/CellOperations.jl:315

Is it possible to remove the periodic boundary conditions at some point using the PeriodicSystems interface?

My use case is to compute and save to disk the unwrapped particle coordinates (without applying periodic boundary conditions) for computing dynamical quantities such as the mean squared displacement.

Maybe multiple dispatch is being abused for UpdateCellLIst!. The code is getting somewhat hard to follow.

Additionally, the parallel flag and the passing of AuxThreaded became in some sense redundant, since parallel=false and the provision of a AuxThreaded preallocated structure does not make sense.

To be studied.

Hello!

Thank you for the help :)

Something similar to this in a sense:

So that one can visualize the grid used for the neighbour search. I know it is uniform cells and the cell width is 2x cutoff, but could not find the cell center coordinates

Kind regards

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

Hello!

In my simulations I often need to define a "box of interest", i.e. any particle which ends up outside this box will be removed from the simulation. I was wondering if there was a way to specify this using this package?

For example for a dam break simulation:

I would like any particle leaving the red box to leave the simulation.

I could properly do this using a package to find the bounding box, then a manual check each time step if any particles has left the domain.

I thought it might be of use for others as well, so that is why I suggest it :)

Kind regards

Hi,

Do you have support for primitive neighbor list? https://wiki.fysik.dtu.dk/ase/_modules/ase/neighborlist.html

I am also using CellListMap in python, memory copy takes long time from list to numpy and torch, is it possible to return numpy array direclty?

Best

Hi, I am reimplementing something very similar to this https://discourse.julialang.org/t/pairwise-computation-slower-than-python-cython-code-balltree-very-slow/62273/24 for two different particle samples. Below the relevant parts of my implementation

using StaticArrays
using LinearAlgebra
using CellListMap

function coordinate_separation(a, b, box_size)
    delta = abs(a - b)
    return (delta > 0.5*box_size ? delta - box_size : delta)*sign(a-b)
end

function separation_vector(a, b, box_size)
    return SVector{3,Float32}(ntuple(i -> coordinate_separation(a[i],b[i],box_size[i]), 3))
end

function map_function(x, y, i, j, d2, output, vx, vy, bin_edges, box_size)
    s_vector = separation_vector(x, y, box_size)
    norm = sqrt(d2)
    bin_id = searchsortedfirst(bin_edges, norm) - 1
    output[1][bin_id] += (LinearAlgebra.dot(vx - vy, s_vector) / norm)
    output[2][bin_id] += 1
    return output
end

function reduce_hist(hist,hist_threaded)
    hist = hist_threaded[1]
    for i in 2:Threads.nthreads()
     hist[1] .+= hist_threaded[i][1]
     hist[2] .+= hist_threaded[i][2]
    end
    return hist
  end



function pairwise_vel_cellist(sample_1::Vector{SVector{3,Float32}},
                                    vel_1::Vector{SVector{3,Float32}}, 
                                    sample_2::Vector{SVector{3,Float32}},
                                    vel_2::Vector{SVector{3,Float32}},
                                    bin_edges::Vector,
                                    box_size::SVector{3},
                                    )
    max_dist = bin_edges[end]
    n_bins = size(bin_edges)
    box = Box(box_size, max_dist)
    cl = CellList(sample_1, sample_2, box)
    print(cl)
    output = (zeros(n_bins),
            zeros(Int32, n_bins))
    map_pairwise!(
        (x, y, i, j, d2, output) -> map_function(x, y, i, j, d2, output, vel_1, vel_2, bin_edges, box_size),
        output, box, cl,
        reduce = reduce_hist,
        show_progress = true
    )

    return output[2], output[1]
end

When I call the function pairwise_vel_cellist on my data I get the following error

ERROR: LoadError: TaskFailedException

    nested task error: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(72945),), b has dims (Base.OneTo(50383),), mismatch at 1")
    Stacktrace:
     [1] promote_shape
       @ ./indices.jl:178 [inlined]

and a much longer trace. But it seems that the complaint is because one set of particles has one dimension (~50k) and the other has ~70k particles.

Have you encountered this issue? I suspect I'm overlooking some detail but I'm pretty new to Julia so it's kind of hard to find.

Thanks in advance!

See:

https://github.com/gromacs/gromacs/blob/612a3f43773a0adbc9e86ab1acef48e5f470beb7/src/gromacs/pbcutil/pbc.cpp#L128

Thank you for creating the high-level PeriodicSystems interface, it has simplified my simulation code.

My simulations typically start without particles and more get added at random times. This makes it tricky to use with PeriodicSystems because currently, it errors if I make a system without particles, see below.

Currently, I am constructing the PeriodicSystem with a few junk particles, and separately keeping track of if it is valid or not. It would simplify my code if PeriodicSystem could properly handle the case of no particles.

using StaticArrays
using CellListMap.PeriodicSystems
Vec2D = SVector{2,Float64}
system = PeriodicSystem(;
    positions=Vec2D[],
    cutoff=0.1,
    unitcell=[1.0, 1.0],
    output=0.0,
)

ERROR: DivideError: integer division error
Stacktrace:
 [1] rem
   @ ./int.jl:289 [inlined]
 [2] set_idxs!(idxs::Vector{UnitRange{Int64}}, n_particles::Int64, nbatches::Int64)
   @ CellListMap ~/.julia/packages/CellListMap/A9riF/src/CellLists.jl:414
 [3] CellListMap.AuxThreaded(cl::CellListMap.CellList{2, Float64}; particles_per_batch::Int64)
   @ CellListMap ~/.julia/packages/CellListMap/A9riF/src/CellLists.jl:392
 [4] AuxThreaded
   @ ~/.julia/packages/CellListMap/A9riF/src/CellLists.jl:373 [inlined]
 [5] #UpdateCellList!#52
   @ ~/.julia/packages/CellListMap/A9riF/src/CellLists.jl:687 [inlined]
 [6] CellListMap.CellList(x::Vector{SVector{2, Float64}}, box::CellListMap.Box{CellListMap.OrthorhombicCell, 2, Float64, Float64, 4}; parallel::Bool, nbatches::Tuple{Int64, Int64})
   @ CellListMap ~/.julia/packages/CellListMap/A9riF/src/CellLists.jl:506
 [7] PeriodicSystem(; positions::Vector{SVector{2, Float64}}, xpositions::Nothing, ypositions::Nothing, unitcell::Vector{Float64}, cutoff::Float64, output::Float64, output_name::Symbol, parallel::Bool, nbatches::Tuple{Int64, Int64}, lcell::Int64)
   @ CellListMap.PeriodicSystems ~/.julia/packages/CellListMap/A9riF/src/PeriodicSystems.jl:165
 [8] top-level scope
   @ REPL[10]:1

Hi! First of all, thanks for putting up together this code. It's been pretty useful for my calculations of galaxy two-point clustering.

I have run a few tests where I compute the PDF of galaxy number densities (which is pretty similar to a histogram of particle distances), and I've found that the number of pairs that are returned by map_pairwise change depending on how many threads are used. I only get sensible results when using a single thread. I have setup up a test repository to exemplify the problem:

https://github.com/epaillas/julia-sandbox

I call Julia from Python by running call_julia.py, which will generate a plot of the galaxy density PDF. I get different results if I change the number of threads in Line 6.

I am using an M1 Macbook Pro 2020 with Big Sur, Python 3.8.8, Julia 1.7.0, and the latest version of CellListMap.

Is there any obvious problem in my code that would cause this discrepancy, or could it be something else?

Hello, thanks a lot for the package. It is really fantastic.

I am currently trying to implement a basic Monte Carlo simulation engine here (https://github.com/edwinb-ai/Metropolis.jl), and I am using CellListMap.jl for computing the interactions between particles.

However, I am having some trouble when defining a type that contains a field that is of type Box. The PR in which I am working is here.
Basically, the problem is that I want to define the following type

mutable struct System{B,T,VT,I}
    xpos::VT    # position of particles
    density::T   # density of the system
    temperature::T    # temperature
    box::B    # simulation box, actually want it to be of type `Box`
    rng::Random.AbstractRNG
    npart::I    # total number of particles
end

I then have a constructor that automatically creates the simulation box, which is this one

function System(
    density::T, temp::T, particles::I, cutoff::T; dims=3, random_init=true, lcell=2
) where {T<:Real,I<:Int}
    box_size = cbrt(T(particles) / density)
    # box = create_box(box_size, dims, cutoff; lcell=lcell)
    box = CellListMap.Box(fill(box_size, dims), cutoff; lcell=lcell)
    rng = Xorshifts.Xoroshiro128Plus()
    xpos = initialize_positions(box_size, rng, particles; random_init=random_init)
    syst = System(xpos, density, temp, box, rng, particles)

    return syst
end

I am using Cthulhu.jl to check for type instabilities, and the one that I get is that the type of box is not concrete, and so there are many allocations happening.
I was wondering if you could point me in the right direction of how I should be declaring the types of Box inside other types, such as the one I described above.

The Box constructor is type-unstable, because of the two types of unit cells and the two dimensions that are possible.

This type instability is benign, as it does not propagate to the performance-critical functions, but anyway it is annoying sometimes, if it is being used from within a larger code base, because it causes some minor allocations.

Document the how to use strip_value.

I found that after updating the InPlaceNeighborList() for serval steps, the result of neighbor!(system) can be different from that of
neighbor(x_new)

using CellListMap, Test
x = rand(SVector{3,Float64}, 10^3);
system = InPlaceNeighborList(x=x, cutoff=0.1, unitcell=[1,1,1], parallel=false)

list = neighborlist!(system)

x_new = rand(SVector{3,Float64}, 10^3);
for _=1:1000
    x_new = rand(SVector{3,Float64}, 10^3);
    update!(system, x_new)
end

@test neighborlist!(system) == neighborlist(x_new, 0.1; unitcell = [1, 1, 1])

I think that may because the change in size of the neighborlist.

we strongly encourage the use of the PeriodicSystems interface

From the docs, it seems that the high-level PeriodicSystems interface is supposed to be used for most purposes.
However, it neither supports non-periodic boxes nor coordinate matrices.

The underlying low-level interface supports both of these, so this is just due to too strict type requirements in the PeriodicSystems interface:

• unitcell must be an abstract vector or matrix and can't be a Limits, even though it's supported by Box.
• xpositions must be a vector of vectors, even though the CellList supports matrices as well.