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
end
end
function applypbc(x::Vector, L)
map!(x->singlepbc(x,L), x,x)
end
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
applypbc(x_new,L)
#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
end
#=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
end
function lennard_jones_force(eps, sigma, r)
force = 4*eps/sigma*(-(sigma/r)^13 + (sigma/r)^7 )
return force
end
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
end
end
end
return a
end
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)
push!(xs,x)
push!(vs,v)
push!(as,a)
end
return xs,vs,as
end
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
...
Stacktrace:
[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.
).