Skip to content

Commit

Permalink
Merge pull request #5 from m3g/sorted
Browse files Browse the repository at this point in the history
merege all changes of sorted branch
  • Loading branch information
lmiq authored Aug 3, 2021
2 parents 91c1e3d + 59121f2 commit 648ec61
Show file tree
Hide file tree
Showing 13 changed files with 1,704 additions and 1,094 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "CellListMap"
uuid = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
authors = ["Leandro Martinez <[email protected]> and contributors"]
version = "0.5.1"
version = "0.5.2"

[deps]
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
Expand Down
50 changes: 19 additions & 31 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -72,18 +72,18 @@ The full code of the examples described here is available at the [examples.jl](h

### Mean difference of coordinates

Computing the mean difference in `x` position between random particles. The closure is used to remove the indexes and the distance of the atoms from the parameters of the input function, as they are not needed in this case.
Computing the mean difference in `x` position between random particles. The closure is used to remove the indexes and the distance of the particles from the parameters of the input function, as they are not needed in this case.

```julia
using CellListMap

# System properties
n = 100_000
N = 100_000
sides = [250,250,250]
cutoff = 10

# Particle positions
x = [ box.sides .* rand(SVector{3,Float64}) for i in 1:n ]
x = [ sides .* rand(3) for i in 1:N ]

# Initialize linked lists and box structures
box = Box(sides,cutoff)
Expand All @@ -94,18 +94,17 @@ f(x,y,sum_dx) = sum_dx + abs(x[1] - y[1])
normalization = N / (N*(N-1)/2) # (number of particles) / (number of pairs)

# Run calculation (0.0 is the initial value)
avg_dx = normalization * map_parwise(
(x,y,i,j,d2,sum_dx) -> (x,y,sum_dx), 0.0, box, cl
avg_dx = normalization * map_pairwise!(
(x,y,i,j,d2,sum_dx) -> f(x,y,sum_dx), 0.0, box, cl
)

```

The example above can be run with `CellListMap.test1()`.

### Histogram of distances

Computing the histogram of the distances between particles (considering the same particles as in the above example). Again,
we use a closure to remove the indexes of the atoms from the computation, because they are not needed. The distance, on the other side, is needed in this example:
we use a closure to remove the indexes of the particles from the computation, because they are not needed. The distance, on the other side, is needed in this example:

```julia
# Function that accumulates the histogram of distances
Expand Down Expand Up @@ -184,7 +183,7 @@ The example above can be run with `CellListMap.test4()`.

### Nearest neighbor

Here we compute the indexes of the atoms that satisfy the minimum distance between two sets of points, using the linked lists. The distance and the indexes are stored in a tuple, and a reducing method has to be defined for that tuple to run the calculation. The function does not need the coordinates of the points, only their distance and indexes.
Here we compute the indexes of the particles that satisfy the minimum distance between two sets of points, using the linked lists. The distance and the indexes are stored in a tuple, and a reducing method has to be defined for that tuple to run the calculation. The function does not need the coordinates of the points, only their distance and indexes.

```julia
# Number of particles, sides and cutoff
Expand Down Expand Up @@ -229,7 +228,7 @@ The example above can be run with `CellListMap.test5()`. The example `CellListMa

### 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.
In this example we compute the complete neighbor list, of all pairs of particles which are closer than the desired cutoff. The implementation returns a vector of tuples, in which each tuple contains the indexes of the particles 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.

```julia
# Function to be evalulated for each pair: push pair if d<cutoff
Expand Down Expand Up @@ -266,9 +265,9 @@ The full example can be run with `CellListMap.test7()`.

## Periodic boundary conditions

**WARNING: The development of general periodic boundary conditions is still premilimary. Non-orthorhombic periodicity is not extensively tested. Next release (0.6) will contain more stable and tested implementation of this feature.**
Triclinic periodic boundary conditions of any kind can be used. However, the input has some limitations for the moment. The lattice vectors must have strictly positive coordinates, and the smallest distance within the cell cannot be smaller than twice the size of the cutoff. An error will be produced if the cell does not satisfy these 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.
Let us illustrate building a two-dimensional cell, 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`):

Expand All @@ -280,11 +279,9 @@ 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)
julia> CellListMap.draw_computing_cell(x,box)
```

<img src=./src/assets/lattice.png>
Expand All @@ -294,28 +291,19 @@ 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 rectangular box (or orthorhombic box, in three-dimensions), with boundaries that exceed the actual system size. This improves the performance of the pairwise computations by avoding the necessity of wrapping coordinates on the main loop (this is an implementation detail only). The resulting box can be visualized with:
109 cells with real particles.
2041 particles in computing box, including images.

```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.
Upon construction of the cell lists, the particles are replicated to fill a rectangular box (or orthorhombic box, in three-dimensions), with boundaries that exceed the actual system size. This improves the performance of the pairwise computations by avoding the necessity of wrapping coordinates on the main loop (this is an implementation detail only).

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> box = Box([ 50. 0. 00.
0. 30. 30.
0. 00. 50. ], 2.)

julia> x = 100*rand(SVector{3,Float64},10000);

Expand Down Expand Up @@ -425,7 +413,7 @@ box = Box(x,box,lcell=2)
cl = CellList(x,box)
map_pairwise!(...)
```
This parameter determines how fine is the mesh of cells. There is a trade-off between the number of cells and the number of particles per cell. For low-density systems, greater meshes are better, because each cell will have only a few atoms and the computations loop over a samller number of cells. For dense systems, it is better to run over more cells with less atoms per cell. It is a good idea to test different values of `lcell` to check which is the optimal choice for your system. Usually the best value is between `lcell=1` and `lcell=6`, but for large and dense systems a larger value may be optimal. For molecular systems with normal densities `lcell=2` is likely the optimal choice. The peformance can be tested using the progress meter, as explained below.
This parameter determines how fine is the mesh of cells. There is a trade-off between the number of cells and the number of particles per cell. For low-density systems, greater meshes are better, because each cell will have only a few particles and the computations loop over a samller number of cells. For dense systems, it is better to run over more cells with less particles per cell. It is a good idea to test different values of `lcell` to check which is the optimal choice for your system. Usually the best value is between `lcell=1` and `lcell=6`, but for large and dense systems a larger value may be optimal. For molecular systems with normal densities `lcell=1` is likely the optimal choice. The peformance can be tested using the progress meter, as explained below.

### Output progress

Expand All @@ -442,7 +430,7 @@ Thus, besides being useful for following the progress of a long run, it is usefu

## Some benchmarks

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:
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. 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.5.1.png>
Expand Down
Loading

2 comments on commit 648ec61

@lmiq
Copy link
Member Author

@lmiq lmiq commented on 648ec61 Aug 3, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/42115

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.2 -m "<description of version>" 648ec6135b36e4b20882f7de4332b3626fd21788
git push origin v0.5.2

Please sign in to comment.