MethodError in 1D Lennard-Jones fluid simulation

Hello people,
I’m super new to Julia (*) and as a starter I am simulating a 1D Lennard-Jones fluid. To do so I’m integrating Newton’s equations with periodic boundary conditions using second order velocity verlet. I tested the integrator with harmonic oscillators and it seems fine. Here’s the code:

using Random
using Plots

function singlepbc(x,L)
    if 0<=x<L
        return x
    elseif x>=L
        return x%L
    elseif x<0
        return x%L + L

function applypbc(x::Vector, L)
    map!(x->singlepbc(x,L), x,x)

function verlet_step(x,v,a,dt,L)
    #new positions 
    N =length(x)
    x_new = zeros(N)    
    a_new = zeros(N)
    v_new = zeros(N)
    x_new .= x .+ v .* dt + (0.5*dt^2)  .* a
    #new accelerations
    a_new = accelerations(x,L)
    #new velocities
    v_new .= v .+ 0.5 .*(a+a_new) .* dt
    return x_new,v_new, a_new 
#=extracts the time serie (eg x(t), v(t),a(t)) from the table (eg xs,vs,as) 
of the particle identified by index =#
function timeserie(index::Int, table)
    timeserie = [elem[index] for elem in table]
    return timeserie
function lennard_jones_force(eps, sigma, r)
    force = 4*eps/sigma*(-(sigma/r)^13 + (sigma/r)^7 )
    return force
function accelerations(x,L)
    #calculates accelerations for 1D lennard-jones particles
    k = 5.0
    m = 1.0
    N = length(x)
    a = zeros(N)
    eps = 1.;
    sigma =1.;

    for i in 1:N
        for j in i+1:N
            rij = x[j]-x[i]
            rij -= L *round(rij/L) #minimum image convention
            if -2.5<rij<2.5
                aij = lennard_jones_force(EPSILON, SIGMA, rij) / m 
                a[i] += aij
                a[j] -= aij #Newton's third law
    return a

function simulate(x0, v0, num_steps, dt, L)
    x = copy(x0)
    v = copy(v0)
    a = accelerations(x,L)
    xs = [x,]
    vs = [v,]
    as = [a,]
    for step in 1:num_steps
        x, v, a = verlet_step( x, v, a, dt, L)
    return xs,vs,as
const global N_PART = 10
const global L_BOX  = 40.
const global EPSILON = 1.
const global SIGMA = 1.

x0 = L_BOX .* rand(N_PART)
v0 = -5. .+ 10*rand(N_PART)
num_steps = 500
dt = 0.002
xs,vs,as = @time simulate(x0,v0,num_steps,dt,L_BOX)
ts = [t for t in 0:dt:dt*num_steps]
particle_trajectory = timeserie(1,xs)
@time plot(ts, [timeserie(k,xs) for k in 1:N_PART]) 

This code seems to work when few collisions occur (e.g. for big L_BOX, or small N_PART). However it usally gives the following error:

ERROR: MethodError: Cannot `convert` an object of type Nothing to an object of type Float64

Closest candidates are:
  convert(::Type{T}, ::T) where T<:Number
   @ Base number.jl:6
  convert(::Type{T}, ::Number) where T<:Number
   @ Base number.jl:7
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number
   @ Base twiceprecision.jl:273

 [1] setindex!(A::Vector{Float64}, x::Nothing, i1::Int64)
   @ Base ./array.jl:969
 [2] map!(f::var"#7#8"{Int64}, dest::Vector{Float64}, A::Vector{Float64})
   @ Base ./abstractarray.jl:3255
 [3] applypbc
   @ ~/Projects/learn_julia/abinfo/lj_temp.jl:15 [inlined]
 [4] verlet_step(x::Vector{Float64}, v::Vector{Float64}, a::Vector{Float64}, dt::Float64, L::Int64)
   @ Main ~/Projects/learn_julia/abinfo/lj_temp.jl:26
 [5] simulate(x0::Vector{Float64}, v0::Vector{Float64}, num_steps::Int64, dt::Float64, L::Int64)
   @ Main ~/Projects/learn_julia/abinfo/lj_temp.jl:77
 [6] top-level scope
   @ ./timing.jl:273 [inlined]
 [7] top-level scope
   @ ~/Projects/learn_julia/abinfo/lj_temp.jl:0

Honestly i don’t understand the error, but my guess is that something blows up when two particles get too close (due to the r^(-13) repulsive force) and a Nothing type value is returned. The fact that the system is in 1D may be relevant to the problem.
It would be great if you could explain me the reason of this error and (if possible) how to fix it, thank you very much!

(*) I couldn’t choose between this section and the “New to Julia” one, I’m sorry if this question is too basic for this section.

The error comes from singlepbc(...) returning nothing because none of the if conditions are satisfied. This happens when x is NaN floating point.

With the first mystery unravelled, now why can x be a NaN is a new question. This probably happens when calculating the Lennard-Jones potential in case r is close to 0.0.

Indeed in the case of small r, the force can blow up. A minimum cap on how close r can get, might solve this problem. In the “real world” particles cannot be closer than some minimum distance (unless fusion is part of the simulation :wink:).

Also, velocities should also be wrapped at L/dt as moving around the box in a time step is unnoticeable to the simulation (as positions are wrapped too).

Thank you for the quick reply. What you say is indeed sensible, I was missing the Nothing returned in case x couldn’t satisfy any condition in applypbc(). I solved the problem by setting force() to a costant below an arbitrary value of r. It’s not the most elegant solution but I guess it’s physically reasonable.

function lennard_jones_force(eps, sigma, r)
    q = sigma/r
    if q < 10^3
        force = 4*eps/sigma*(-(q)^13 + (q)^7 )
        force = -4*eps/sigma*10^39
    return force

EDIT: I didn’t see you suggested the same solution to the problem. Thank you very much for the help, this solution seems to work. What you say about the velocities is interesting, I am almost fully convinced: the only doubt i have is that it seems I’m throwing away kinetic energy.