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
# 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