# Monte Carlo particle simulations are very slow

Hello everyone,

I’m still learning Julia by doing some standard physics simulations in Julia. I used to have all of my code base in C and Fortran, and just recently attempted to migrate everything to Julia.

In particular, I’m doing a very simple Monte Carlo simulation of a Lennard-Jones fluid in the NVT ensemble. The thing is that in C this code would be very fast, finishing in almost 2-3 minutes, at most. But using the `ProgressMeter.jl` I’m finding out that the same exact code will take up to 1 hour in Julia.

I’m trying to figure out what is wrong with my code, where the bottleneck is, but I’m stuck, any help would be greatly appreciated. Here is the full code of my implementation. For now, I’m interested in solving the performance issues.

``````using Pkg
Pkg.activate("./")
Pkg.instantiate()

module MonteCarlo

using StaticArrays
using Random
using RandomNumbers.Xorshifts
using ProgressMeter

mutable struct System{V <: Real}
N::Integer
L::V
ρ::V
press::V
β::V
δr::V
rc2::V
corr::V
rng
ener::V
end

function energy(pos::AbstractArray, syst::System)
total_energy = 0
x = view(pos, :, 1)
y = view(pos, :, 2)
z = view(pos, :, 3)

@inbounds for i = 1:(syst.N - 1)
@fastmath for j = (i + 1):syst.N
xij = x[i] - x[j]
yij = y[i] - y[j]
zij = z[i] - z[j]

# Periodic boundaries
xij -= syst.L * round(xij / syst.L)
yij -= syst.L * round(yij / syst.L)
zij -= syst.L * round(zij / syst.L)

# Compute distance
Δpos2 = xij*xij + yij*yij + zij*zij

if Δpos2 < syst.rc2
r6 = 1.0 / Δpos2^3
# Evaluate the Lennard-Jones potential
total_energy += 4.0 * (r6^2 - r6)
end
end
end
end

function init!(pos::AbstractArray, syst::System)
n3 = 2
ix = 0
iy = 0
iz = 0
x = view(pos, :, 1)
y = view(pos, :, 2)
z = view(pos, :, 3)

# Find the lowest perfect cube, n3, greater than or equal to the
# number of particles
while n3^3 < syst.N
n3 += 1
end

@inbounds for i = 1:syst.N
x[i] = (ix + 0.5) * syst.L / n3
y[i] = (iy + 0.5) * syst.L / n3
z[i] = (iz + 0.5) * syst.L / n3
ix += 1

if ix == n3
ix = 0
iy += 1
if iy == n3
iy = 0
iz += 1
end
end
end
end

function metropolis!(pos::AbstractArray, syst::System)
# Create a new possible displacement
randvec = @SVector rand(3)
δpos = @. syst.δr * (0.5 - randvec)

# Choose a random particle
rng_part = rand(syst.rng, 1:syst.N)
# Save the previous energy
posold = pos[rng_part, :]

# Move that particle
posnew = posold .+ δpos
pos[rng_part, :] = posnew

# Apply periodic boundary conditions
@. posnew -= syst.L * round(posnew / syst.L)

# Compute the energy of the new configuration
newener = energy(pos, syst)

# Always accept lower energy configurations
if newener < syst.ener
syst.ener = newener
return 1
# If not possible, apply the Metropolis acceptance criteria
elseif rand(syst.rng) < exp(syst.β * (syst.ener - newener))
syst.ener = newener
return 1
else
pos[rng_part, :] = posold
return 0
end
end

function run()
N = 500 # number of particles
ρ = 0.86 # reduced density
volume = N / ρ # reduced volume
Lbox = ∛volume # Cubic box length
T = 0.85 # Reduced temperature
P = 1.266 # Reduced pressure
rc = 3.0 # Cut radius
# Long range corrections for the Lennard-Jones potential
corr = 1.0 / rc^3
corr = (corr^3 / 3.0) - corr
corr *= 8.0 * π * ρ / 3.0
# Create a rng
rng = Xorshift1024Plus(123456)
# Initialize the system object
syst = System(N, Lbox, ρ, P, 1 / T, 0.4, rc^2, corr, rng, 0.0)

# Create the positions vector
positions = zeros(Float64, syst.N, 3)
# Initialize the positions as grid
init!(positions, syst)
# Compute the initial energy of the system
# Should be ≈ -5.4518
syst.ener = energy(positions, syst)
println(syst.ener / syst.N)

# Run values
cycles = 100000
attempts = 0
accepted = 0
total_energy = 0
ratio = 0

# Equilibration steps
@showprogress for i = 1:50000
attempts += 1
mc = metropolis!(positions, syst)
accepted += mc
ratio = accepted / attempts
# Save that attempt's energy
total_energy += syst.ener
end
println("Equilibration energy per particle: \$(total_energy / N)")
end

# Run the full simulation
run()

end

``````

One problem is that `N` is an `Integer` (not an `Int`). The former is an abstract type that will thus cause type instability and performance issues. Also, making `Pos` be an array of arrays will be more efficient since you then don’t need to use weird views. Also, what’s `Lcaja`? you use it but never define it. Is it supposed to be `Lbox` (I’m guessing you translated some of your code from spanish?)?

2 Likes

In addition to the problems that @Oscar_Smith mentioned, in your `energy` function you initialise `total_energy` as an `Int`, and inside the loop it becomes a `Float64`. In the end the compiler can’t know exactly which type `total_energy` has, which leads to poor performance.

In this specific case you could initialise `total_energy = 0.0` to fix the issue.

You may want to take a look at the performance tips in the Julia docs, and in particular to this section which describes this precise form of type instability.

2 Likes

This used to be a problem, but the compiler is usually good at figuring this out nowadays. Nonetheless, it’s good practice to initialize variables to the correct type, for performance and to work well with exotic types.

Yes, sorry for that, I’m not a native English language speaker and all of my code was in Spanish. I just fixed that, thanks for noticing.

I see, how can I do that? Could it be

``````positions = Array{Vector{Float64},N}(undef, 3)
``````

which would initialize 3 `Vector{Float64}` arrays of length `N` with `N` being the total number of particles?

1 Like

You want

``````positions = Vector{Vector{Float64}}(undef, 3)
``````

as you allocated a `N` dimensional array.

Actually you should really use an array of SVector from the StaticArrays package for coordinate vectors. Vastly more efficient for fixed size (3-component) arrays.

7 Likes

`System.rng` is also untyped.

Thank you everyone, I have achieved the performance I was hoping for. I removed the `views` and started using `StaticArrays`, I have also typed all that was untyped and it has helped greatly.

I have used most of the suggestions here, so I want to thank everyone who replied, this community is fantastic.

5 Likes

That’s great to hear!

But, in case it was unclear:

The typing that affects performance is on the struct fields. Explicit typing of function arguments does not affect performance.

5 Likes