Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
@vchuravy

This comment has been minimized.

Copy link
Owner Author

@vchuravy vchuravy commented Nov 30, 2014

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment