Lennard-Jones_MD

function ljconfig(N, per=0.2)
    # mesh-grid-like construction
    X = (0:N) * ones(1,N)
    Y = X'
    # write all particles into a single array
    x = [X[:]', Y[:]']
    x = x + per * rand(size(x))
end

N=50

x = ljconfig(N)

function lj!(dJ :: Vector{Float64}, r :: Matrix{Float64}, k :: Int64 ,n :: Int64)
  for i in 1:length(dJ)
    @inbounds	dJ[i] = r[i, k] - r[i, n]
  end
  
  s = sumabs2(dJ)

  J = s^(-6) - 2 * s^(-3) # 47%
  Δ = -12.0 * (s^(-7) - s^(-4)) # 47%
  for i in 1:length(dJ)
    @inbounds dJ[i] = Δ * dJ[i]
  end
  return J
end

function lj_rewrite_math!(dJ :: Vector{Float64}, r :: Matrix{Float64}, k :: Int64 ,n :: Int64)
  for i in 1:length(dJ) # 6%
    @inbounds	dJ[i] = r[i, k] - r[i, n] # 9%
  end

  s = sumabs2(dJ) # 6%
  t = 1.0 /s^3 # 19%
  t² = t^2 # 10%
  J = t² - 2.0 * t # 4%
  Δ = -12.0 * (t² - t) / s # 35%
  for i in 1:length(dJ)
    @inbounds dJ[i] = Δ * dJ[i] # 7%
  end
  return J
end

function lj(r)
    s = sumabs2(r)
    J = s^(-6) - 2 * s^(-3)
    dJ = -12. * (s^(-7) - s^(-4)) * r
    return J, dJ
end

function lj_pretty(x)
  E = 0.0
  dE = zeros(size(x))
  for n = 1:(size(x,2)-1), k = (n+1):size(x,2)
    r = x[:,k]-x[:,n] # 31%
    J, dJ = lj(r) # 29%
    E += J
    dE[:, k] += dJ # 20%
    dE[:, n] -= dJ # 19%
  end
  return E, dE
end

function lj_pretty_1(x)
  E = 0.0
  dE = zeros(size(x))
  for n = 1:(size(x,2)-1), k = (n+1):size(x,2)
    r = x[:,k]-x[:,n] # 50%
    J, dJ = lj(r) # 50%
    E += J
    for i in 1:size(dE, 1)
      @inbounds begin
        Δ = dJ[i]
        dE[i, k] += Δ
        dE[i, n] -= Δ
      end
    end
  end
  return E, dE
end

function lj_pretty_2(x)
  E = 0.0
  dE = zeros(size(x))
  dJ = zeros(size(x, 1))
  for n = 1:(size(x,2)-1), k = (n+1):size(x,2)
    E += lj!(dJ, x, k, n) # 100%
    for i in 1:size(dE, 1)
      @inbounds begin
        Δ = dJ[i]
        dE[i, k] += Δ
        dE[i, n] -= Δ
      end
    end
  end
  return E, dE
end

function lj_pretty_2_rewrite_math(x)
  E = 0.0
  dE = zeros(size(x))
  dJ = zeros(size(x, 1))
  for n = 1:(size(x,2)-1), k = (n+1):size(x,2)
    E += lj_rewrite_math!(dJ, x, k, n) # 100%
    for i in 1:size(dE, 1)
      @inbounds begin
        Δ = dJ[i]
        dE[i, k] += Δ
        dE[i, n] -= Δ
      end
    end
  end
  return E, dE
end

E, dE = lj_pretty(x)
E1, dE1 = lj_pretty_1(x)
E2, dE2 = lj_pretty_2(x)
E3, dE3 = lj_pretty_2_rewrite_math(x)

@time lj_pretty(x)
@time lj_pretty_1(x)
@time lj_pretty_2(x)
@time lj_pretty_2_rewrite_math(x)

E == E1 == E2 # true
isapprox(E, E3) # true
all(dE .== dE1) # true
all(dE .== dE2) # true
all(map(isapprox, dE, dE3)) # true

Timings:

elapsed time: 6.109372319 seconds (2834019224 bytes allocated, 30.61% gc time)
elapsed time: 3.367907747 seconds (1378030424 bytes allocated, 26.37% gc time)
elapsed time: 1.051502035 seconds (41088 bytes allocated)
elapsed time: 0.109610573 seconds (41088 bytes allocated)