With 15 lines I could write a code to perform a particle simulation with periodic boundary conditions, a langevin thermostat, and a quadratic potential between the particles, and produce an animation:

```
using Plots ; ENV["GKSwstype"]="nul"
const N, τ, Δt, λ, T, k = 100, 1000, 0.01, 1e-3, 0.26, 1e-6
const x, v, f = -0.5 .+ rand(3,N), -0.01 .+ 0.02*randn(3,N), zeros(3,N)
wrap(x,y) = (x-y) > 0.5 ? (x-y)-1 : ( (x-y) < -0.5 ? (x-y)+1 : (x-y) )
anim = @animate for t in 1:τ
f .= 0
for i in 1:N-1, j in i+1:N
f[:,i] .+= wrap.(x[:,i],x[:,j]) .- λ .* v[:,i]
f[:,j] .+= wrap.(x[:,j],x[:,i]) .- λ .* v[:,j]
end
x .= wrap.(x .+ v*Δt .+ (f/2)*Δt^2,zeros(3))
v .= v .+ f*Δt .+ sqrt.(2*λ*k*T*Δt)*randn()
scatter(x[1,:],x[2,:],x[3,:],label="",lims=(-0.5,0.5),aspect_ratio=1,framestyle=:box)
end
gif(anim,"anim.gif",fps=10)
```

As Tamas mentioned, if we go to a few more lines, this could include an actual Lennard-Jones potential and be parallelized. And that without using any package besides Plots.