Accelerate pairwise Lennard-Jones force computation

General Usage Performance
82

EDIT: I got the concept to work and can do the calculation as I want in theory :slight_smile: I noticed though a few things:

  1. In the NeighbourList or the mapping procedure, the particle it self will not be added to the neighbour list. This means you will never see (i,i,d2), which in most cases is fine, but for my calculation I for example need the property of i too.

  2. It seems like that for some reason (maybe my wrong use of box, using “limits(x)” now), that not all neighbouring points are found - this is based on my comparison with the NearestNeighbours package which seems to find all neighbours as I would expect

EDIT: I just realized that my second point is properly because that the CellListMap is smart and only saves one entry for distance! So if i and j are neighbours with d2 = 4, then it will only save (i,j,4) and not (j,i,4) which makes sense!

Just Complaining

I am really struggling to figure out how to do the calculation as you mentioned:

function f(x,y,i,j,d2,output)
    output += # your function of coordinates x and y, indexes i, and j, and distance sqrt(d2)
    return output
end

I don’t know from where I have to provide these things? Do I need to do the neighbour list first, then I would have x, y, i, j and d2?

How would I go about doing it as you described without outputting neighbour list?

All the examples also mention box and celllist, do I need to generate this? Can it not be done automatically based on the points I input or?

Sorry for the many questions, I am just quite stuck.

Kind regards

1 Like
83

Indeed, it does not, but then the distance is zero and nothing has to be computed, a straight loop over the particles should solve that.

I will write a more detailed step by step and add it to the manual. For now the documentation is relying heavily on the examples.

84

Yep exactly, that is also how I ended up doing it

And sure, would help a lot of people out - the forces example made me finally grasp it

1 Like
85

Let me know, in private if you prefer, of any difficulties, feedbacks, etc. If you find it slow or something else. All feedback is appreciated.

1 Like
86

Do you also intend on providing functionality with your package to analyse triplet functions for example, or more complex particle properties that depend on pair distance vectors such as (averaged) Steinhardt order parameters, or performing a Voronoi tessellation?

87

The function that is mapped receives the particle relative positions for each pair, so something as you suggest may be possible to implement as an application of it. I am not familiar with the details of those, so I would need to study them to be completely sure.

(Distribution fictions are certainly an application, I already use it for that in the ComplexMixtures.jl package)

1 Like
88

Ah that makes sense, I agree that the latter two might indeed be easy to implement using your package. I might have a go at that at some point.

I took a look at your package CompleMixtures.jl. Cool stuff! Pair distribution functions are indeed a very natural application of this. I am looking forward to seeing more of your work.

2 Likes
89

Nice discussion! I just found that sorting the pairlist could improve the performance mainly because the r = x - y in the first line in the function force_pair. After I commented out this line and replace r with a fixed random static vector, sorting the pairlist contributes little to the performance. Is there any possible explantation?

using BenchmarkTools
using StaticArrays
using LoopVectorization

@inline @fastmath function force_pair!(x,y,i,j,d2,ε,σ6,f,r)
    # r = y - x
    d8 = (d2^2)^2
    σ6d8 = σ6*inv(d8)
    d23inv = inv(d2^3)
    dudr = -12ε * σ6d8 * (σ6*d23inv - 1)*r
    @inbounds f[i] += dudr
    @inbounds f[j] -= dudr
    return f
end

function compute_force!(f,x,y,pairlist,distances,r)
    σ = 3.28018
    ε = 0.0441795
    σ6 = σ^6
    for ipair in 1:length(pairlist)
        i = pairlist[ipair][1]
        j = pairlist[ipair][2]
        d2 = distances[ipair]
        f = force_pair!(x[i],y[j],i,j,d2,ε,σ6,f,r)
    end
    return f
end

# data
n = 10_000 # number of particles
npairs = 8026946 # number of pairs within cutoff
pairlist = [ SVector{2,Int}(rand(1:n),rand(1:n)) for _ in 1:npairs ]
sort!(pairlist)
distances = [ 144*rand() for _ in 1:npairs ]
x = [ 12*rand(SVector{2,Float64}) for _ in 1:n ]
y = [ 12*rand(SVector{2,Float64}) for _ in 1:n ]

# Initialize force
f = zeros(SVector{2,Float64},length(x))
r = rand(SVector{2,Float64})

# Benchmark
@benchmark compute_force!($f,$x,$y,$pairlist,$distances,$r)

The preformance is 27ms. After I commented out the sort!(pairlist), the performance is 27ms, too.