# Monte Carlo particle simulations are very slow

**URL:** <https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239>\
**Category:** Performance\
**Created:** [May 10, 2020, 6:45pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239 "2020-05-10T18:45:23Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [May 10, 2020, 6:45pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/1 "2020-05-10T18:45:23Z")

</div>

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`](https://github.com/timholy/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.

```julia
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
    # Add the tail corrections
    return total_energy + (syst.N * syst.corr)
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

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 10, 2020, 6:54pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/2 "2020-05-10T18:54:23Z")

</div>

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?)?

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [May 10, 2020, 7:35pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/3 "2020-05-10T19:35:15Z")

</div>

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](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-changing-the-type-of-a-variable-1) which describes this precise form of type instability.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 10, 2020, 7:47pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/4 "2020-05-10T19:47:33Z")

</div>

> [@jipolanco](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [May 10, 2020, 8:03pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/5 "2020-05-10T20:03:32Z")

</div>

> [@Oscar\_Smith](#):
>
> 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?)

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.

> [@Oscar\_Smith](#):
>
> making `Pos` be an array of arrays will be more efficient since you then don’t need to use weird views

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 10, 2020, 8:04pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/6 "2020-05-10T20:04:49Z")

</div>

You want

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

```

as you allocated a `N` dimensional array.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 10, 2020, 8:32pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/7 "2020-05-10T20:32:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 10, 2020, 9:35pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/8 "2020-05-10T21:35:33Z")

</div>

> [@Oscar\_Smith](#):
>
> One problem is that `N` is an `Integer`

`System.rng` is also untyped.

---

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [May 11, 2020, 4:54pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/9 "2020-05-11T16:54:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 11, 2020, 5:15pm UTC](https://discourse.julialang.org/t/monte-carlo-particle-simulations-are-very-slow/39239/10 "2020-05-11T17:15:40Z")

</div>

That’s great to hear!

But, in case it was unclear:

> [@edwinb-ai](#):
>
> I have also typed all that was untyped and it has helped greatly.

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