# Solid boundary for ABM model

Hello everyone,

I’m building a model where particles roam around randomly in 3d space. At the moment issue is that I want particles to be contained in space though what I see is that it escapes boundary and comes on the other side eg: CellListMap.jl · Agents.jl Here you can see some particles after reaching boundary line comes from other side. What I want is they bounce back inside after coming in contact with compartment boundary rather going to other side. I am using a elastic collision function but the issue is that it seems to work on particles but I don’t know how it can be also set for boundary of the container.

Did you set periodic = false?

1 Like

it helps with restricting particles inside container but at the moment they are sticking to boundary. Is there any option to bounce them back? or it will happen automatically once elastic collision is applied to particles?

You need to create walls that change the particle velocities, and implement that explicitly, computing that for each particle.

Not using periodic boundary conditions, in CellListMap, doesn’t have the desired effect, it makes the system literary without boundaries.

Can you suggest (link or docs or some example already implement this in Julia ) to read to understand implementation? because Im not sure if Agents.jl has any in build functions for this .

code to make reflection on a wall (here called “obstacle”) with a given normal vector n and particle velocity p.vel. In the case of a wall n is a constant, it doesn’t depend on the position p.pos.

3 Likes

Here is a modification of the example with two options to set walls, rigid or “soft”. Rigid walls invert the velocity when the particle reaches the walls. “Soft” walls apply a force that is dependent on the penetration of the particle into the wall, and can be made stiffer for increased strength of such a force. Both “work” reasonably. Soft bouncing tends to be numerically more stable.

The particle positions are initialized within the walls, which are contained in the “periodic” system, which is large enough for the particles not interact with images. CellListMap only handles non-periodicity in this way.

Here is the code:

Agents.jl + CellListMap.jl
# # Integrating Agents.jl with CellListMap.jl

# @raw html
# <video width="auto" controls autoplay loop>
# <source src="../celllistmap.mp4" type="video/mp4">
# </video>
# 

# This example illustrates how to integrate Agents.jl with
# [CellListMap.jl](https:://github.com/m3g/CellListMap.jl), to accelerate the
# computation of short-ranged (within a cutoff) interactions in 2D and 3D continuous
# spaces. CellListMap.jl is a package that allows the computation of pairwise interactions
# using an efficient and parallel implementation of [cell lists](https://en.wikipedia.org/wiki/Cell_lists).

# ## The system simulated
#
# The example will illustrate how to simulate a set of particles in 2 dimensions, interacting
# through a simple repulsive potential of the form:
#
# $U(r) = k_i k_j\left[r^2 - (r_i+r_j)^2\right]^2~~~\textrm{for}~~~r \leq (r_i+r_j)$
#
# $U(r) = 0.0~~~\textrm{for}~~~r \gt (r_i+r_j)$
#
# where $r_i$ and $r_j$ are the radii of the two particles involved, and
# $k_i$ and $k_j$ are constants associated to each particle. The potential
# energy function is a smoothly decaying potential with a maximum when
# the particles overlap.
#
# Thus, if the maximum sum of radii between particles is much smaller than the size
# of the system, cell lists can greatly accelerate the computation of the pairwise forces.
#
# Each particle will have different radii and different repulsion force constants and masses.
using Agents

# Below we define the Particle type, which represents the agents of the
# simulation. The Particle type, for the ContinuousAgent{2,Float64} space, will have additionally
# an id and pos (position) and vel (velocity) fields, which are automatically added
# by the @agent macro.
@agent struct Particle(ContinuousAgent{2,Float64})
k::Float64 # repulsion force constant
mass::Float64
end
PropParticle(; vel, r, k, mass) = (vel, r, k, mass)

# ## Required and data structures for CellListMap.jl
#
# We will use the high-level ParticleSystem interface (requires version ≥0.9.0):
using CellListMap

# Two auxiliary arrays will be created on model initialization, to be passed to
# the ParticleSystem data structure:
#
# 1. positions: CellListMap requires a vector of (preferentially) static vectors as the positions
#    of the particles. To avoid creating this array on every call, a buffer to
#    which the agent.pos positions will be copied is stored in this data structure.
# 2. forces: In this example, the property to be computed using CellListMap.jl is
#    the forces between particles, which are stored here in a Vector{<:SVector}, of
#    the same type as the positions. These forces will be updated by the map_pairwise!
#    function.
#
# Additionally, the computation with CellListMap.jl requires the definition of a cutoff,
# which will be twice the maximum interacting radii of the particles, and the geometry of the
# the system, given by the unitcell of the periodic box. For non-periodic systems,
# use unitcell = nothing.
#
# More complex output data, variable system geometries and other options are supported,
# according to the [CellListMap](https://m3g.github.io/CellListMap.jl/stable/ParticleSystem/)
# user guide.
#
# ## Model initialization
# We create the model with a keyword-accepting function as is recommended in Agents.jl.
# The keywords here control number of particles and sizes.
function initialize_bouncing(;
number_of_particles=10_000,
sides=SVector(600.0, 600.0),
dt=0.001,
parallel=true
)
## the walls will be at 50.0 and 550.0
## initial particle positions are within the walls
positions = [50.0 .+ 500.0 .* rand(SVector{2,Float64}) for _ in 1:number_of_particles]

## We will use CellListMap to compute forces, with similar structure as the positions
forces = similar(positions)

## Space and agents: <= no need for periodic boundary conditions here,
## although this is sort of a hack in combination with CellListMap.
space2d = ContinuousSpace(sides; periodic=false)

## Initialize CellListMap particle system
system = ParticleSystem(
positions=positions,
unitcell=sides,
output=forces,
output_name=:forces, # allows the system.forces alias for clarity
parallel=parallel,
)

## define the model properties
## The system field contains the data required for CellListMap.jl
properties = (
dt=dt,
number_of_particles=number_of_particles,
system=system,
)
model = StandardABM(Particle,
space2d;
agent_step!,
model_step!,
agents_first = false,
properties=properties
)

## Create active agents
for id in 1:number_of_particles
pos = positions[id]
prop_particle = PropParticle(
r = (0.5 + 0.9 * rand()) * max_radius,
k = 10 + 20 * rand(), # random force constants
mass = 10.0 + 100 * rand(), # random masses
vel = 100 * randn(SVector{2}) # initial velocities)
)
end

return model
end

# ## Computing the pairwise particle forces
# To follow the CellListMap interface, we first need a function that
# computes the force between a single pair of particles. This function
# the periodic boundary conditions), x and y, their indices in the
# array of positions, i and j, the squared distance between them, d2,
# the forces array to be updated and the model properties.
#
# Given these input parameters, the function obtains the properties of
# each particle from the model, and computes the force between the particles
# as (minus) the gradient of the potential energy function defined above.
#
# The function *must* return the forces array, to follow the CellListMap API.
#
function calc_forces!(x, y, i, j, d2, forces, model)
p_i = model[i]
p_j = model[j]
d = sqrt(d2)
if d ≤ (p_i.r + p_j.r)
dr = y - x # x and y are minimum-image relative coordinates
fij = 2 * (p_i.k * p_j.k) * (d2 - (p_i.r + p_j.r)^2) * (dr / d)
forces[i] += fij
forces[j] -= fij
end
return forces
end

# The model_step! function will use CellListMap to update the
# forces for all particles. The first argument of the call is
# the function to be computed for each pair of particles, which closes-over
# the model data to call the calc_forces! function defined above.
#
function model_step!(model::ABM)
## Update the pairwise forces at this step
map_pairwise!(
(x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, forces, model),
model.system,
)
return nothing
end

# ## Update agent positions and velocities
# The agent_step! function will update the particle positions and velocities,
# given the forces, which are computed in the model_step! function. A simple
# Euler step is used here for simplicity. Finally, the positions stored in the ParticleSystem
# structure are updated.
function agent_step!(agent, model::ABM)
id = agent.id
dt = abmproperties(model).dt
## Retrieve the forces on agent id
f = model.system.forces[id]
# Choose wall type
walltype = :rigid
##
## Add bouncing on walls force F(x) = -k * Δx with large k
##
if walltype == :soft
fx, fy = 0.0, 0.0
k = 10^6
if agent.pos[1] ≤ 50.0
fx +=  k * (50.0 - agent.pos[1])
end
if agent.pos[1] ≥ 550.0
fx += k * (550.0 - agent.pos[1])
end
if agent.pos[2] ≤ 50.0
fy += k * (50.0 - agent.pos[2])
end
if agent.pos[2] ≥ 550.0
fy += k * (550.0 - agent.pos[2])
end
# sum wall forces to the particle forces
f += SVector(fx, fy)
end
a = f / agent.mass
## Update positions and velocities
v = agent.vel + a * dt
x = agent.pos + v * dt + (a / 2) * dt^2
##
## rigid walls
##
xx, xy = x
vx, vy = v
if walltype == :rigid
xx, xy = x
vx, vy = v
if xx ≤ 50.0
xx = 50.0
vx < 0.0 && (vx = -vx)
end
if xx ≥ 550.0
xx = 550.0
vx > 0.0 && (vx = -vx)
end
if xy ≤ 50.0
xy = 50.0
vy < 0.0 && (vy = -vy)
end
if xy ≥ 550.0
xy = 550.0
vy > 0.0 && (vy = -vy)
end
end
x = SVector(xx, xy)
agent.vel = SVector(vx, vy)
move_agent!(agent, x, model)
## !!! IMPORTANT: Update positions in the ParticleSystem
model.system.positions[id] = agent.pos
return nothing
end

# ## The simulation
# Finally, the function below runs an example simulation, for 1000 steps.
function simulate(model=nothing; nsteps=1_000, number_of_particles=10_000)
if isnothing(model)
model = initialize_bouncing(number_of_particles=number_of_particles)
end
Agents.step!(model, nsteps)
end
# Which should be quite fast
#model = initialize_bouncing()
#@time simulate(model)

# and let's make a nice video with less particles,
# to see them bouncing around. The marker size is set by the
# radius of each particle, and the marker color by the
# corresponding repulsion constant.

using CairoMakie
#CairoMakie.activate!() # hide
model = initialize_bouncing(number_of_particles=1000)
abmvideo(
"celllistmap.mp4", model;
framerate=20, frames=200, dt=5,
title="Softly bouncing particles with CellListMap.jl",
agent_size=p -> p.r,
agent_color=p -> p.k
)

# @raw html
# <video width="auto" controls autoplay loop>
# <source src="../celllistmap.mp4" type="video/mp4">
# </video>
# 



3 Likes

I will try to implement this modification on my model. I don’t think there should be too much issue with one more dimension.

Thank you

1 Like

I think it is a good idea to add a new function in Agents.jl which updates the position and velocity of the agent accordingly, shouldn’t be too hard @Sahil_Khan if you want to contribute that feature given that you are working on that in your model

I’ve edited my reply: i have given incorrect code by mistake. The correct reflection code is actually only 2 lines of code.

1 Like

code:

# # Integrating Agents.jl with CellListMap.jl

# @raw html
# <video width="auto" controls autoplay loop>
# <source src="../celllistmap.mp4" type="video/mp4">
# </video>
# 

# This example illustrates how to integrate Agents.jl with
# [CellListMap.jl](https:://github.com/m3g/CellListMap.jl), to accelerate the
# computation of short-ranged (within a cutoff) interactions in 2D and 3D continuous
# spaces. CellListMap.jl is a package that allows the computation of pairwise interactions
# using an efficient and parallel implementation of [cell lists](https://en.wikipedia.org/wiki/Cell_lists).

# ## The system simulated
#
# The example will illustrate how to simulate a set of particles in 2 dimensions, interacting
# through a simple repulsive potential of the form:
#
# $U(r) = k_i k_j\left[r^2 - (r_i+r_j)^2\right]^2~~~\textrm{for}~~~r \leq (r_i+r_j)$
#
# $U(r) = 0.0~~~\textrm{for}~~~r \gt (r_i+r_j)$
#
# where $r_i$ and $r_j$ are the radii of the two particles involved, and
# $k_i$ and $k_j$ are constants associated to each particle. The potential
# energy function is a smoothly decaying potential with a maximum when
# the particles overlap.
#
# Thus, if the maximum sum of radii between particles is much smaller than the size
# of the system, cell lists can greatly accelerate the computation of the pairwise forces.
#
# Each particle will have different radii and different repulsion force constants and masses.
using Agents

# Below we define the Particle type, which represents the agents of the
# simulation. The Particle type, for the ContinuousAgent{2,Float64} space, will have additionally
# an id and pos (position) and vel (velocity) fields, which are automatically added
# by the @agent macro.
@agent struct Particle(ContinuousAgent{3,Float64})
k::Float64 # repulsion force constant
mass::Float64
end
PropParticle(; vel, r, k, mass) = (vel, r, k, mass)

# ## Required and data structures for CellListMap.jl
#
# We will use the high-level ParticleSystem interface (requires version ≥0.9.0):
using CellListMap

# Two auxiliary arrays will be created on model initialization, to be passed to
# the ParticleSystem data structure:
#
# 1. positions: CellListMap requires a vector of (preferentially) static vectors as the positions
#    of the particles. To avoid creating this array on every call, a buffer to
#    which the agent.pos positions will be copied is stored in this data structure.
# 2. forces: In this example, the property to be computed using CellListMap.jl is
#    the forces between particles, which are stored here in a Vector{<:SVector}, of
#    the same type as the positions. These forces will be updated by the map_pairwise!
#    function.
#
# Additionally, the computation with CellListMap.jl requires the definition of a cutoff,
# which will be twice the maximum interacting radii of the particles, and the geometry of the
# the system, given by the unitcell of the periodic box. For non-periodic systems,
# use unitcell = nothing.
#
# More complex output data, variable system geometries and other options are supported,
# according to the [CellListMap](https://m3g.github.io/CellListMap.jl/stable/ParticleSystem/)
# user guide.
#
# ## Model initialization
# We create the model with a keyword-accepting function as is recommended in Agents.jl.
# The keywords here control number of particles and sizes.
function initialize_bouncing(;
number_of_particles=10_000,
sides=SVector(600.0, 600.0, 600.0),
dt=0.001,
parallel=true
)
## the walls will be at 50.0 and 550.0
## initial particle positions are within the walls
positions = [50.0 .+ 500.0 .* rand(SVector{3,Float64}) for _ in 1:number_of_particles]

## We will use CellListMap to compute forces, with similar structure as the positions
forces = similar(positions)

## Space and agents: <= no need for periodic boundary conditions here,
## although this is sort of a hack in combination with CellListMap.
space2d = ContinuousSpace(sides; periodic=false)

## Initialize CellListMap particle system
system = ParticleSystem(
positions=positions,
unitcell=sides,
output=forces,
output_name=:forces, # allows the system.forces alias for clarity
parallel=parallel,
)

## define the model properties
## The system field contains the data required for CellListMap.jl
properties = (
dt=dt,
number_of_particles=number_of_particles,
system=system,
)
model = StandardABM(Particle,
space2d;
agent_step!,
model_step!,
agents_first = false,
properties=properties
)

## Create active agents
for id in 1:number_of_particles
pos = positions[id]
prop_particle = PropParticle(
r = (0.5 + 0.9 * rand()) * max_radius,
k = 10 + 20 * rand(), # random force constants
mass = 10.0 + 100 * rand(), # random masses
vel = 100 * randn(SVector{3}) # initial velocities)
)
end

return model
end

# ## Computing the pairwise particle forces
# To follow the CellListMap interface, we first need a function that
# computes the force between a single pair of particles. This function
# the periodic boundary conditions), x and y, their indices in the
# array of positions, i and j, the squared distance between them, d2,
# the forces array to be updated and the model properties.
#
# Given these input parameters, the function obtains the properties of
# each particle from the model, and computes the force between the particles
# as (minus) the gradient of the potential energy function defined above.
#
# The function *must* return the forces array, to follow the CellListMap API.
#
function calc_forces!(x, y, i, j, d2, forces, model)
p_i = model[i]
p_j = model[j]
d = sqrt(d2)
if d ≤ (p_i.r + p_j.r)
dr = y - x # x and y are minimum-image relative coordinates
fij = 2 * (p_i.k * p_j.k) * (d2 - (p_i.r + p_j.r)^2) * (dr / d)
forces[i] += fij
forces[j] -= fij
end
return forces
end

# The model_step! function will use CellListMap to update the
# forces for all particles. The first argument of the call is
# the function to be computed for each pair of particles, which closes-over
# the model data to call the calc_forces! function defined above.
#
function model_step!(model::ABM)
## Update the pairwise forces at this step
map_pairwise!(
(x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, forces, model),
model.system,
)
return nothing
end

# ## Update agent positions and velocities
# The agent_step! function will update the particle positions and velocities,
# given the forces, which are computed in the model_step! function. A simple
# Euler step is used here for simplicity. Finally, the positions stored in the ParticleSystem
# structure are updated.
function agent_step!(agent, model::ABM)
id = agent.id
dt = abmproperties(model).dt
## Retrieve the forces on agent id
f = model.system.forces[id]
# Choose wall type
walltype = :rigid
##
## Add bouncing on walls force F(x) = -k * Δx with large k
##
if walltype == :soft
fx, fy, fz = 0.0, 0.0, 0.0
k = 10^6
if agent.pos[1] ≤ 50.0
fx +=  k * (50.0 - agent.pos[1])
end
if agent.pos[1] ≥ 550.0
fx += k * (550.0 - agent.pos[1])
end
if agent.pos[2] ≤ 50.0
fy += k * (50.0 - agent.pos[2])
end
if agent.pos[2] ≥ 550.0
fy += k * (550.0 - agent.pos[2])
end
if agent.pos[3] ≤ 50.0
fz += k * (50.0 - agent.pos[3])
end
if agent.pos[3] ≥ 550.0
fz += k * (550.0 - agent.pos[3])
end
# sum wall forces to the particle forces
f += SVector(fx, fy, fz)
end
a = f / agent.mass
## Update positions and velocities
v = agent.vel + a * dt
x = agent.pos + v * dt + (a / 2) * dt^2
##
## rigid walls
##
xx, xy, xz = x
vx, vy, xz = v
if walltype == :rigid
xx, xy, xz = x
vx, vy, vz = v
if xx ≤ 50.0
xx = 50.0
vx < 0.0 && (vx = -vx)
end
if xx ≥ 550.0
xx = 550.0
vx > 0.0 && (vx = -vx)
end
if xy ≤ 50.0
xy = 50.0
vy < 0.0 && (vy = -vy)
end
if xy ≥ 550.0
xy = 550.0
vy > 0.0 && (vy = -vy)
end
if xz ≤ 50.0
xz = 50.0
vz < 0.0 && (vz = -vz)
end
if xz ≥ 550.0
xz = 550.0
vz > 0.0 && (vz = -vz)
end
end
x = SVector(xx, xy, xz)
agent.vel = SVector(vx, vy, vz)
move_agent!(agent, x, model)
## !!! IMPORTANT: Update positions in the ParticleSystem
model.system.positions[id] = agent.pos
return nothing
end

# ## The simulation
# Finally, the function below runs an example simulation, for 1000 steps.
function simulate(model=nothing; nsteps=1_000, number_of_particles=10_000)
if isnothing(model)
model = initialize_bouncing(number_of_particles=number_of_particles)
end
Agents.step!(model, nsteps)
end
# Which should be quite fast
#model = initialize_bouncing()
#@time simulate(model)

# and let's make a nice video with less particles,
# to see them bouncing around. The marker size is set by the
# radius of each particle, and the marker color by the
# corresponding repulsion constant.

using CairoMakie
#CairoMakie.activate!() # hide
model = initialize_bouncing(number_of_particles=100)
abmvideo(
"celllistmap.mp4", model;
framerate = 20, frames = 50, dt=5,
title="Softly bouncing particles with CellListMap.jl",
agent_size=p -> p.r,
agent_color=p -> p.k
)

# @raw html
# <video width="auto" controls autoplay loop>
# <source src="../celllistmap.mp4" type="video/mp4">
# </video>
# 


Video output:

1 Like

sorry if its naive but i didn’t get input type for obstacle. For particle it will be agent which can be extracted from model I believe but i’m not sure I understand obstacle. I mean like what input it will take?

Hey, I had some questions related to these two functions as I want to modify them to also include attraction after collision:

function calc_forces!(x, y, i, j, d2, forces, model)
p_i = model[i]
p_j = model[j]
d = sqrt(d2)
if d ≤ (p_i.r + p_j.r)
dr = y - x # x and y are minimum-image relative coordinates
fij = 2 * (p_i.k * p_j.k) * (d2 - (p_i.r + p_j.r)^2) * (dr / d)
forces[i] += fij
forces[j] -= fij
end
return forces
end

# The model_step! function will use CellListMap to update the
# forces for all particles. The first argument of the call is
# the function to be computed for each pair of particles, which closes-over
# the model data to call the calc_forces! function defined above.
#
function model_step!(model::ABM)
## Update the pairwise forces at this step
map_pairwise!(
(x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, forces, model),
model.system,
)
return nothing
end

1. here, x and y are as same as p_i.pos, p_j.pos ?
2. From where d2 comes in function ?
3. As in model_step! function we are not giving any input. how does it picks put x,y, d2 etc ? from model element?
4. Can you also provide function name used for force calculation. so i can check out on web?

Thank you

They are, but wrapped one to each other according to periodic boundary conditions (if these are used).

d2 is the squared distance between the particles (sqrt(d2) is the distance between them). It is provided by the interface of CellListMap because it is computed already internally, so you can skip computing it again.

The data is all inside the model.system object. What the map_pairwise! is saying is: For all pairs of particles that are closer to each other than the cutoff, apply the function calc_forces!. The interface is read like:

(x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, forces, model),


on the left side are the default parameters of the interface: the function to be applied for each pair can receive the coordinates, x and y of the particles, their indices i, j, the squared distance between them d2, and an output variable, here named forces. In this case, this function will be mapped to the call to calc_forces! because this one receives an extra parameter, model.

That specific function does not have a “name”. It is the local repulsion function used in Packmol. It is a soft and differentiable repulsion function. But for any specific problem you probably want something tailored for your problem.

ps: you might find this notebook interesting.

Thank you:). So, if I understand correctly x, y are coordinates for particles. So, in case of 3d it should be x, y, z ? and same changes should be made for force calculation?

Do you think it will be much easier to implement just elastic collision? because in my case i want particles to stop based on type of particles when they interact with each other and stay stick.

Oh, no, x and y are the vectors of particle coordinates of each particle. Thus, each is a two- or three-dimensional vector.

That force function, specifically, is generic, it works for two or three dimensions, because dr = y - x is a vector, thus fij is a vector, and forces[i] += fij is adding vector fij to the force on atom i, where forces is a vector of 2D- or 3D vectors.

like so:

julia> using StaticArrays

julia> x = SVector(1.0, 2.0, 3.0); y = SVector(4.0, 5.0, 6.0);

julia> dr = x - y
3-element SVector{3, Float64} with indices SOneTo(3):
-3.0
-3.0
-3.0

1 Like

okay…make sense why it worked with 3d implementation without any changes

1 Like