Merge pull request #4 from m3g/pbc

merged pdc

Project.toml

@@ -1,10 +1,11 @@
 name = "CellListMap"
 uuid = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 authors = ["Leandro Martinez <[email protected]> and contributors"]
-version = "0.5.0"
+version = "0.5.1"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -15,6 +16,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 DocStringExtensions = "0.8"
 Parameters = "0.12"
 ProgressMeter = "1.6"
+Setfield = "0.7"
 StaticArrays = "1"
 julia = "1"
 

README.md

@@ -1,6 +1,6 @@
 # CellListMap.jl
 
-This package is for computing short-ranged particle interactions or any other property that is dependent on the distances between pairs of three-dimensional particles, within a cutoff. It maps a function to be computed pairwise using cell lists using, for the moment, orthorhombic periodic boundary conditions. Parallel and serial implementations can be used. 
+This package is for computing short-ranged particle interactions or any other property that is dependent on the distances between pairs of three-dimensional particles, within a cutoff. It maps a function to be computed pairwise using cell lists using periodic boundary conditions of any type. Parallel and serial implementations can be used. 
 
 It allows the fast computation of any quantity from the pairs that are within the desired cutoff, for example an average distance or an histogram of distances, forces, potentials, minimum distances, etc., as the examples below illustrate. This is done by passing the function to be evaluated as a parameter of the `map_pairwise!` function. 
 
@@ -14,7 +14,8 @@ It allows the fast computation of any quantity from the pairs that are within th
       3. [Gravitational potential](#gravitational-potential) 
       4. [Gravitational force](#gravitational-force) 
       5. [Nearest neighbor](#nearest-neighbor) 
-      6. [Neighbor list](#neighbor-list) 
+      6. [Neighbour list](#neighbor-list) 
+      7. [Periodic boundary conditions](#periodic-boundary-conditions)
 4. [Parallelization splitting and reduction](#parallelization-splitting-and-reduction)
       1. [Custom reduction functions](#custom-reduction-functions) 
 5. [Preallocating auxiliary arrays: threaded output and cell lists](#preallocating-auxiliary-arrays-threaded-output-and-cell-lists)
@@ -226,7 +227,7 @@ mind = map_pairwise!(
 
 The example above can be run with `CellListMap.test5()`. The example `CellListMap.test6()` of [examples.jl](https://github.com/m3g/CellListMap.jl/blob/8661ae692abf3f44094f1fc41076c464300729b6/src/examples.jl#L219) describes a similar problem but *without* periodic boundary conditions. Depending on the distribution of points it is a faster method than usual ball-tree methods. 
 
-### Neighbor list
+### Neighbour list
 
 In this example we compute the complete neighbor list, of all pairs of atoms which are closer than the desired cutoff. The implementation returns a vector of tuples, in which each tuple contains the indexes of the atoms and the corresponding distance. The empty `pairs` output array will be split in one vector for each thread, and reduced with a custom reduction function. 
 
@@ -263,6 +264,67 @@ pairs = map_pairwise!(
 
 The full example can be run with `CellListMap.test7()`. 
 
+## Periodic boundary conditions
+
+Periodic boundary conditions of any kind can be used. Let us illustrate its use with a two-dimensional case, for easier visualization. A matrix of column-wise lattice vectors is provided in the construction of the box, and that is all. 
+
+Here, the lattice vectors are `[1,0]` and `[0.5,1]` (and we illustrate with `cutoff=0.1`): 
+
+```julia
+julia> box = Box([ 1.0  0.5
+                     0  1.0 ], 0.1);
+
+julia> x = 10*rand(SVector{2,Float64},1000);
+```
+We have created random coordinates for `1000` particles, that are not necessarily wrapped according to the periodic boundary conditions. We can see the coordinates in the minimum image cell with:
+```julia
+julia> p = [ CellListMap.wrap_to_first(x,box) for x in x ];
+
+julia> using Plots
+
+julia> scatter(Tuple.(p),aspect_ratio=1,framestyle=:box,label=:none)
+```
+
+<img src=./src/assets/lattice.png>
+
+The construction of the cell list is, as always, done with:
+
+```julia
+julia> cl = CellList(x,box)
+CellList{2, Float64}
+  90 cells with real particles.
+  2065 particles in computing box, including images.
+
+```
+
+Upon construction of the cell lists, the particles are replicated to fill a square box (or cubic box, in three-dimenions). This improves the performance of the pairwise computations by avoding the necessity of wrapping coordinates on the main loop. The resulting box can be visualized with:
+
+```julia
+julia> p = CellListMap.view_celllist_particles(cl,box);
+
+julia> scatter(Tuple.(p),aspect_ratio=1,framestyle=:box,label=:none)
+```
+<img src=./src/assets/replicated.png>
+
+These are the particles that will be considered for the calculations. 
+
+In summary, to use arbitrary periodic boundary conditions, just initialize the box with the matrix of lattice vectors. In three dimensions, for example, one could use:
+
+```julia
+julia> box = Box([ 50.  0. 20. 
+                    0.  0. 30.          
+                    0. 30. 50. ], 10.)
+
+julia> x = 100*rand(SVector{3,Float64},10000);
+
+julia> p = [ CellListMap.wrap_to_first(x,box) for x in x ];
+
+julia> scatter(Tuple.(p),aspect_ratio=1,framestyle=:box,label=:none)
+```
+to work with an arbitrary 3D lattice, Which in this case looks like:
+
+<img src=./src/assets/3Dlattice.png>
+
 ## Parallelization splitting and reduction
 
 The parallel execution requires the splitting of the computation among threads, obviously. Thus, the output variable must be split and then reduced to avoid concurrency. To control these steps, set manually the `output_threaded` and `reduce` optional input parameters of the `map_pairwise!` function. 
@@ -381,9 +443,9 @@ Thus, besides being useful for following the progress of a long run, it is usefu
 The goal here is to provide a good implementation of cell lists. We compare it with the implementation of the nice cython/python [halotools](https://github.com/astropy/halotools) package, in the computation of an histogram of mean pairwise velocities. These tests are implemented in the [halotools.jl](https://github.com/m3g/CellListMap.jl/blob/main/src/halotools.jl) file. Currently, the `CellListMap.jl` is as fast for dense systems, and scales linearly and parallelizes well for increasing number of particles, with constant density:
 
 
-<img src=https://github.com/lmiq/PairVelocities/blob/main/data/cd_v0.4.8.png>
+<img src=https://github.com/lmiq/PairVelocities/blob/main/data/cd_v0.5.1.png>
 
-<img src=https://github.com/lmiq/PairVelocities/blob/main/data/cv_v0.4.8.png>
+<img src=https://github.com/lmiq/PairVelocities/blob/main/data/cv_v0.5.1.png>
 
 The full test is available [at this](https://github.com/lmiq/PairVelocities) repository. And we kindly thank [Carolina Cuesta](https://github.com/florpi) for providing the example. 
 

src/CellListMap.jl

@@ -468,7 +468,7 @@ function ranges_of_replicas(cutoff,nc,unit_cell,unit_cell_max::SVector{2,T}) whe
   )
   r_min, r_max = _ranges_of_replicas(
     nc,
-    SVector{3,Int}(-1,-1),
+    SVector{2,Int}(-1,-1),
     unit_cell,
     cell_vertices
   )
@@ -554,6 +554,17 @@ first unit cell with all-positive coordinates.
   p = wrap_cell_fraction(x,cell)
   return cell*p
 end
+"""
+
+```
+wrap_to_first(x::SVector{N,T},box::Box)
+```
+
+Wraps the coordinates of point `x` such that the returning coordinates are in the
+first unit cell with all-positive coordinates, given the `Box` structure.
+
+"""
+@inline wrap_to_first(x,box::Box) = wrap_to_first(x,box.unit_cell)
 
 """
 
@@ -633,23 +644,52 @@ end
 """
 
 ```
-neighbour_cells(lcell::Int)
+neighbour_cells(box::Box{N,T,M}) where {N,M}
 ```
 
-Function that returns the iterator of the cartesian indices of the cell that must be 
-evaluated if the cells have sides of length `cutoff/lcell`. 
+Function that returns the iterator of the cartesian indices of the cells that must be 
+evaluated (forward, i. e. to avoid repeated interactions) 
+if the cells have sides of length `box.cutoff/box.lcell`. `N` can be
+`2` or `3`, for two- or three-dimensional systems.
 
 """
-function neighbour_cells(lcell)
+function neighbour_cells(box::Box{3,T,9}) where T
+  @unpack lcell = box
   nb = Iterators.flatten((
     CartesianIndices((1:lcell,-lcell:lcell,-lcell:lcell)),
     CartesianIndices((0:0,1:lcell,-lcell:lcell)),
     CartesianIndices((0:0,0:0,1:lcell))
   ))
   return nb
 end
-neighbour_cells_all(lcell) = 
-  CartesianIndices((-lcell:lcell,-lcell:lcell,-lcell:lcell))
+function neighbour_cells(box::Box{2,T,4}) where T
+  @unpack lcell = box
+  nb = Iterators.flatten((
+    CartesianIndices((1:lcell,-lcell:lcell)),
+    CartesianIndices((0:0,1:lcell))
+  ))
+  return nb
+end
+
+"""
+
+```
+neighbour_cells_all(box::Box{N,T,M}) where {N,M}
+```
+
+Function that returns the iterator of the cartesian indices of all neighbouring
+cells of a cell if the cells have sides of `box.cutoff/box.lcell`. `N` can be
+`2` or `3`, for two- or three-dimensional systems.
+
+"""
+function neighbour_cells_all(box::Box{3,T,9}) where T  
+  @unpack lcell = box
+  return CartesianIndices((-lcell:lcell,-lcell:lcell,-lcell:lcell))
+end
+function neighbour_cells_all(box::Box{2,T,4}) where T  
+  @unpack lcell = box
+  return CartesianIndices((-lcell:lcell,-lcell:lcell))
+end
 
 """
 
@@ -869,7 +909,7 @@ function inner_loop!(f,box,icell,cl::CellList,output)
     i = pᵢ.index
   end
 
-  for jcell in neighbour_cells(box.lcell)
+  for jcell in neighbour_cells(box)
     output = cell_output!(f,box,cell,cl,output,cell.cartesian+jcell)
   end
 
@@ -912,7 +952,6 @@ end
 # Serial version for cross-interaction computations
 #
 function map_pairwise_serial!(f::F, output, box::Box, cl::CellListPair; 
-  parallel::Bool=false,
   show_progress=show_progress
 ) where {F}
   show_progress && (p = Progress(length(cl.small),dt=1))
@@ -929,7 +968,6 @@ end
 function map_pairwise_parallel!(f::F1, output, box::Box, cl::CellListPair;
   output_threaded=output_threaded,
   reduce::F2=reduce,
-  parallel::Bool=false,
   show_progress=show_progress
 ) where {F1,F2}
   show_progress && (p = Progress(length(cl.small),dt=1))
@@ -947,10 +985,10 @@ end
 # Inner loop of cross-interaction computations
 #
 function inner_loop!(f,output,i,box,cl::CellListPair)
-  @unpack unit_cell, nc, lcell, cutoff_sq = box
+  @unpack unit_cell, nc, cutoff_sq = box
   xpᵢ = wrap_to_first(cl.small[i],unit_cell)
   ic = particle_cell(xpᵢ,box)
-  for neighbour_cell in neighbour_cells_all(lcell)
+  for neighbour_cell in neighbour_cells_all(box)
     jc = cell_linear_index(nc,neighbour_cell+ic)
     pⱼ = cl.large.fp[jc]
     j = pⱼ.index

src/testing.jl

@@ -60,20 +60,21 @@ julia> scatter(x,y,label=nothing,xlims=(-10,180),ylims=(-10,180))
 ```
 
 """
-function view_celllist_particles(cl::CellList{N,T}) where {N,T}
-  @unpack ncwp, cwp, ncp, fp, np = cl
-  x = Vector{SVector{N,T}}(undef,ncp[1])
+function view_celllist_particles(cl::CellList{N,T},box::Box) where {N,T}
+  @unpack nc = box
+  @unpack cwp, ncp, fp, np = cl
+  x = Vector{SVector{N,T}}(undef,ncp)
   ip = 0
-  for i in 1:ncwp[1]
+  for p in cl.fp
+    p.index == 0 && continue
     ip += 1
-    p = fp[cwp[i].icell]
     x[ip] = p.coordinates
     while np[p.index].index > 0
       ip += 1
       x[ip] = p.coordinates
       p = np[p.index]
     end
   end
-  return ([x[i][j] for i in 1:ncp[1]] for j in 1:N)
+  return [SVector{N,T}(ntuple(j -> x[i][j],N)) for i in 1:ncp]
 end