# Solid boundary for ABM model

**URL:** <https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984>\
**Category:** Modelling & Simulations\
**Tags:** question, agents, modelling, abm, agentsjl\
**Created:** [May 8, 2024, 7:11am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984 "2024-05-08T07:11:33Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 9, 2024, 11:51pm UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/22 "2024-12-09T23:51:12Z")

</div>

Basically you have to replace the computation of forces of the code above by [this one](https://github.com/m3g/2021_FortranCon/blob/46559659ddd45dd17f129ce5bfce930dc8b6abce/celllistmap_vs_namd/simulate.jl#L40).

But the issue is that with a potential energy that becomes very repulsive at short distances you’ll need an [initial set of coordinates without overlaps](https://m3g.github.io/Packmol.jl/stable/#Packmol.pack_monoatomic!) and a short time-step.

Out of curiosity: what’s the motivation for this?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [December 10, 2024, 12:25am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/23 "2024-12-10T00:25:15Z")

</div>

Thanks for the advice!

I am examining how the physical environment affects interaction patterns of social agents. I would like to leverage simple particle physics to allow agents to move through the environment and interact with certain constraints imposed. The model will build upon this example [Continuous space social distancing · Agents.jl Example Zoo](https://juliadynamics.github.io/AgentsExampleZoo.jl/dev/examples/social_distancing/). I thought CellListMap would speed up the simulation, but I am concerned about the impact of particles moving through other particles. Maybe I should stick with the slower setup in the model zoo above.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 10, 2024, 12:34am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/24 "2024-12-10T00:34:41Z")

</div>

Cool. CellListMap has nothing to do with the particles overlapping or not. It replaces the loop over interacting pairs, and the interactions can be of any type. In that example elastic collusions. That can definitely work and be accelerated, if necessary.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [December 10, 2024, 10:35am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/25 "2024-12-10T10:35:03Z")

</div>

@lmiq, thanks for all of your help. I think I am close to getting the non-overlapping example working. I was hoping to get your feedback on my changes, and how to resolve an issue where the force components are NaN when agents reach the boundary. Here is the print out when the error occurs:

```julia
current position [550.0, 550.0] new position [NaN, NaN]
current veolocity [-5.705101313537676e12, -1.00689134264408e12] new velocity [NaN, NaN]
current fource [NaN, NaN]

```

I attached the current code below:

> **Summary**
>
> ````julia
> # # 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.
> ##################################################################################################################################
> cd(@ __DIR__ )
> using Pkg
> Pkg.activate("../../..")
> using Agents
> using CairoMakie
> using CellListMap
> using FastPow
> using Packmol
> using StaticArrays
> 
> # 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})
> r::Float64 # radius
> 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):
> 
> # 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(;
> n_particles = 10_000,
> sides = SVector(600.0, 600.0),
> dt = 0.001,
> max_radius = 10.0,
> ε = 0.0441795,
> σ = 2 * 1.64009,
> parallel = true
> )
> ## the walls will be at 50.0 and 550.0
> ## initial particle positions are within the walls
> positions = [500.0 .* rand(SVector{2, Float64}) for _ = 1:n_particles]
> unitcell = [500, 500]
> pack_monoatomic!(positions, unitcell, max_radius)
> # add [50,50] to make min x = 50, and min y = 50
> map!(x -> x + [50, 50], positions, positions)
> ## 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,
> cutoff = 2 * max_radius,
> 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,
> n_particles,
> system,
> ε,
> σ
> )
> model = StandardABM(Particle,
> space2d;
> agent_step!,
> model_step!,
> agents_first = false,
> properties = properties
> )
> 
> ## Create active agents
> for id = 1:n_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)
> )
> add_agent!(pos, Particle, model, prop_particle...)
> 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
> # receives the positions of the two particles (already considering
> # 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
> 
> function calc_forces!(x, y, i, j, d2, forces, model)
> ε = model.ε
> σ = model.σ
> r = y - x
> @fastpow dudr = -12 * ε * (σ^12 / d2^7 - σ^6 / d2^4) * r
> @inbounds forces[i] = forces[i] + dudr
> @inbounds forces[j] = forces[j] - dudr
> 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)
> if any(isnan, x)
> println("current position $(agent.pos) new position $x")
> println("current veolocity $(agent.vel) new velocity $v")
> println("current force $f")
> end
> 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, n_particles = 10_000)
> if isnothing(model)
> model = initialize_bouncing(n_particles = n_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(n_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>
> # ```
> 
> ````

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 10, 2024, 11:30am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/26 "2024-12-10T11:30:36Z")

</div>

I don´t see anything clearly wrong. Have you tried: 1) using smaller initial velocities? 2) reducing the time-step?

From the output you showed, the issue seems more that the velocity is too high at that step, which is probably a consequence of the simulation becoming unstable because of the issues above.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [December 10, 2024, 11:40am UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/27 "2024-12-10T11:40:51Z")

</div>

Thanks for the suggestion. I don’t have much intuition for the appropriate range of parameter values, especially for ε and σ. So I used the values in your example on GitHub. Reducing the velocities by half seems to have fixed the issue. Interestingly, the particles still pass through each other. Do you have any idea why that might be the case?

> **Summary**
>
> ````julia
> # # 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.
> ##################################################################################################################################
> cd(@ __DIR__ )
> using Pkg
> Pkg.activate("../../..")
> using Agents
> using CairoMakie
> using CellListMap
> using FastPow
> using Packmol
> using StaticArrays
> 
> # 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})
> r::Float64 # radius
> 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):
> 
> # 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(;
> n_particles = 10_000,
> sides = SVector(600.0, 600.0),
> dt = 0.001,
> max_radius = 10.0,
> ε = 0.0441795,
> σ = 2 * 1.64009,
> parallel = true
> )
> ## the walls will be at 50.0 and 550.0
> ## initial particle positions are within the walls
> positions = [500.0 .* rand(SVector{2, Float64}) for _ = 1:n_particles]
> unitcell = [500, 500]
> pack_monoatomic!(positions, unitcell, max_radius)
> # add [50,50] to make min x = 50, and min y = 50
> map!(x -> x + [50, 50], positions, positions)
> ## 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,
> cutoff = 2 * max_radius,
> 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,
> n_particles,
> system,
> ε,
> σ
> )
> model = StandardABM(Particle,
> space2d;
> agent_step!,
> model_step!,
> agents_first = false,
> properties = properties
> )
> 
> ## Create active agents
> for id = 1:n_particles
> pos = positions[id]
> prop_particle = PropParticle(
> #r = (0.5 + 0.9 * rand()) * max_radius,
> r = max_radius,
> k = 10 + 20 * rand(), # random force constants
> mass = 10.0 + 100 * rand(), # random masses
> vel = 50 * randn(SVector{2}) # initial velocities)
> )
> add_agent!(pos, Particle, model, prop_particle...)
> 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
> # receives the positions of the two particles (already considering
> # 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
> 
> function calc_forces!(x, y, i, j, d2, forces, model)
> ε = model.ε
> σ = model.σ
> r = y - x
> @fastpow dudr = -12 * ε * (σ^12 / d2^7 - σ^6 / d2^4) * r
> @inbounds forces[i] = forces[i] + dudr
> @inbounds forces[j] = forces[j] - dudr
> 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)
> if any(isnan, x)
> println("current position $(agent.pos) new position $x")
> println("current veolocity $(agent.vel) new velocity $v")
> println("current force $f")
> end
> 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, n_particles = 10_000)
> if isnothing(model)
> model = initialize_bouncing(n_particles = n_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(n_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>
> # ```
> 
> ````

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 10, 2024, 12:54pm UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/28 "2024-12-10T12:54:37Z")

</div>

I think it might be only an issue of visualization, but the fact is that the code needed some adjustments:

1. Here `sigma` plays the role of the radius (this ia a [lennard-jones](https://en.wikipedia.org/wiki/Lennard-Jones_potential) potential). Thus the `r` and `k` parameters were not really used for the computation of the forces, and could cause confusion.

2. Thus, the number of particles, sigma, and size of the “box” have to talk to each other. The simulation examples I posted above was 3D and had parameters with physical meaning. Here we went 2D and thus to see something reasonable the parameters have to be adjusted.

3. Finally, the `abmvideo` function is too slow to generate a nice video of this type of system. Thus, I print intermediate steps in a trajectory file and created the video with the `VMD` molecular visualization tool. And that for a trajectory of 10^5 steps, which runs in 15 seconds here, generating a more meaningful particle motions.

Thus, I modified a bit the code with those considerations in mind:

> <https://gist.github.com/lmiq/a944b09d43eeedfac05ea1ef56bba898>

And got this video:

[![](https://global.discourse-cdn.com/julialang/original/3X/e/4/e4e7a40d46b895efe14d081ac16728b8ebca513e.jpeg "Lennard-Jones particles in 2D bwith Agents.jl and CellListMap.jl") ](https://www.youtube.com/watch?v=Kw0mSsjVFjo)

It seems everything is as it should be. But we need to be careful with these kind of simulations, in many aspects.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [December 10, 2024, 1:07pm UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/29 "2024-12-10T13:07:42Z")

</div>

Awesome!

I realized I overlooked a critical detail and would like to seek some general advice if that is OK. Where would I include code similar to `transmit!` below?

```julia
function sir_model_step!(model)
    r = model.interaction_radius
    for (a1, a2) in interacting_pairs(model, r, :nearest)
        transmit!(a1, a2, model.reinfection_probability)
        elastic_collision!(a1, a2, :mass)
    end
end

```

Do you recommend adding it before `map_pairwise` and defining some distance threshold that would trigger the function call?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 10, 2024, 1:18pm UTC](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984/30 "2024-12-10T13:18:07Z")

</div>

That should be part of what we are calling here `calc_forces!`. The `calc_forces!` function is the function called for each pair of particles closer to each other than a cutoff, and thus there is where one should deal with changing properties that depend on the proximity of the particles, for performance.

But note that there are again other things to take into consideration: `interacting_pairs(model, r; :nearest)` I suppose is iterating over the **closest** particle of each other particle. That’s one specific type of iteration. What we do in the example above is iterating over all particles that are within a cutoff radius.

If the collision is strictly elastic, there are no forces. Then, what we call `calc_forces!` here should instead change the velocities of the particles according to the elastic collision (elastic-collision simulations are tricky, note the remark in the Agents.jl example as well).

[Previous page](https://discourse.julialang.org/t/solid-boundary-for-abm-model/113984.md?page=1)