Skip to content
Snippets Groups Projects
Commit 2bcc5aa4 authored by Jesse Huhtala's avatar Jesse Huhtala
Browse files

continuing the implementation of energy minimizer

parent 43b362d9
No related branches found
No related tags found
No related merge requests found
...@@ -2,7 +2,8 @@ using Optim ...@@ -2,7 +2,8 @@ using Optim
using Plots using Plots
nel = 10 nel = 10
nfuns = 60 nfuns = 60
r = range(0.01,stop=8.0,length=1000) mu = 2.8
r = range(0.01,stop=7.99,length=1000)
L = 8.0 L = 8.0
bet = zeros(1000) bet = zeros(1000)
rho0 = zeros(1000) rho0 = zeros(1000)
...@@ -77,33 +78,143 @@ function gradient_tf_vw!(G,x) ...@@ -77,33 +78,143 @@ function gradient_tf_vw!(G,x)
G[i] = G[i] -1.0/18.0 * sum(ddbs .* p)*(r[2]-r[1]) G[i] = G[i] -1.0/18.0 * sum(ddbs .* p)*(r[2]-r[1])
end end
end end
function energy_ng(x) function energy_tf_vw_ln(x)
p = psi(x)
ps = dpsi(x)
p2 = ddpsi(x)
en = -0.5*sum(exp.(0.5*p) .* (0.25*ps.^2.0 .* exp.(0.5*p) +0.5*p2 .* exp.(0.5*p) ) )*
(r[2]-r[1])
en = en + sum(0.25*exp.(2*p) - mu*exp.(p))*(r[2]-r[1])
en = en + pi^2.0/6.0 * sum(exp.(3*p))*(r[2]-r[1])
end
function gradient_tf_vw_ln!(G,x)
p = psi(x) p = psi(x)
ps = dpsi(x) ps = dpsi(x)
en = sum(abs.(2 .* p .* ps).^2.0/p.^2.0)*(r[2]-r[1]) p2 = ddpsi(x)
en = en + sum(0.25*p.^4.0+2.0*p.^2.0)*(r[2]-r[1]) for i=1:nfuns
en = en - 2.0*sum(p.^2.0 .* g)*(r[2]-r[1]) bs = basisfunc(i)
nonloc = sum(p.^2.0 .* 2.0 .* (ps) .* p)*(r[2]-r[1])./(p.^2.0) dbs = dbasisfunc(i)
ddbs = ddbasisfunc(i)
#Gradient contribution arising from Thomas-Fermi
G[i] = sum(pi^2.0/2.0 * bs .* exp.(3*p))*(r[2]-r[1])
#Gradient contribution from mu and delta interaction
G[i] = G[i] + sum(0.5*bs .* exp.(2*p)-mu*bs .* exp.(p))
#Gradient contribution arising from von Weizsäcker
G[i] = G[i] - 0.5 * sum(0.5*bs .* ( 0.25*ps.^2.0 .* exp.(0.5*p) +0.5*p2 .* exp.(0.5*p)))*
(r[2]-r[1])
G[i] = G[i] - 0.5 * sum(exp.(0.5*p) .* (0.5*ps .* dbs .* exp.(0.5*p) +
0.25*0.5*ps.^2.0 .* bs .* exp.(0.5*p) + 0.5*ddbs .* exp.(0.5*p) +
0.25*p2 .* bs .* exp.(0.5*p)) )*(r[2]-r[1])
end
end
function energy_ng(x)
p = psi(x)
dp = dpsi(x)
p2 = ddpsi(x)
en = -2.0 * sum(g .* p.^2.0)*(r[2]-r[1])
en = en - sum(mu * p.^2.0)
nonloc = sum(( 2.0*p.^2.0 .* p .* dp))*(r[2]-r[1])
nonloc = nonloc ./ (p.^2.0)
en = en + sum(nonloc)*(r[2]-r[1]) en = en + sum(nonloc)*(r[2]-r[1])
en = en -1.0/2.0 * sum(p .* p2)*(r[2]-r[1])
en = en + sum(0.25 * p.^4.0)*(r[2]-r[1])
return en return en
end end
function gradient_ng!(G,x)
p = psi(x)
p2 = ddpsi(x)
dp = dpsi(x)
for i=1:nfuns
bs = basisfunc(i)
dbs = dbasisfunc(i)
ddbs = ddbasisfunc(i)
#Gradient contribution arising from delta-interaction and mu
G[i] = sum(bs .* p.^3.0 - 2.0*mu* bs .* p)*(r[2]-r[1])
#Gradient contribution arising from von Weizsäcker
G[i] = G[i] -1.0/2.0 * sum(p2 .* bs)*(r[2]-r[1])
G[i] = G[i] -1.0/2.0 * sum(ddbs .* p)*(r[2]-r[1])
#Gradient contribution arising from the non-local part
nonloc1 = (sum(4.0 * p .* bs .* p .* dp)*(r[2]-r[1])) ./(p.^2.0)
nonloc2 = (sum(2.0 * p.^2.0 .* dbs .* p)*(r[2]-r[1])) ./(p.^2.0)
nonloc3 = (sum(2.0 * bs .* p.^2.0 .* dp)*(r[2]-r[1])) ./(p.^2.0)
nonloc4 = (sum(-4.0* p.^2.0 .* p .* dp )*(r[2]-r[1])) ./(p.^3.0)
G[i] = G[i] + sum(nonloc1 .+ nonloc2 .+ nonloc3 .+ nonloc4 .* bs)*(r[2]-r[1])
end
end
function energy_ng_nonl(x)
p = psi(x)
dp = dpsi(x)
p2 = ddpsi(x)
en = -2.0 * sum(g .* p.^2.0)*(r[2]-r[1])
en = en - sum(mu * p.^2.0)
#nonloc = sum(( 2.0*p.^2.0 .* p .* dp))*(r[2]-r[1])
#nonloc = nonloc ./ (p.^2.0)
#en = en + sum(nonloc)*(r[2]-r[1])
en = en -1.0/2.0 * sum(p .* p2)*(r[2]-r[1])
en = en + sum(0.25 * p.^4.0)*(r[2]-r[1])
end
function gradient_ng_nonl!(G,x)
p = psi(x)
p2 = ddpsi(x)
for i=1:nfuns
bs = basisfunc(i)
ddbs = ddbasisfunc(i)
#Gradient contribution arising delta interaction
G[i] = sum(bs .* p.^3.0 - 2*mu*bs .* p)*(r[2]-r[1])
#Gradient contribution arising from von Weizsäcker
G[i] = G[i] -1.0/2.0 * sum(p2 .* bs)*(r[2]-r[1])
G[i] = G[i] -1.0/2.0 * sum(ddbs .* p)*(r[2]-r[1])
#Gradient contribution arising from part of VP
G[i] = G[i] - 4*sum(g .* bs .* p)*(r[2]-r[1])
end
end
function energy_ng_ln(x) function energy_ng_ln(x)
p = psi(x) p = psi(x)
ps = dpsi(x) ps = dpsi(x)
en = sum(abs.(ps .* exp.(p)).^2.0/exp.(p))*(r[2]-r[1]) p2 = ddpsi(x)
en = en+sum(0.25*exp.(p).^2.0-2.0*exp.(p))*(r[2]-r[1]) en = -0.5*sum(exp.(0.5*p) .* (0.25*ps.^2.0 .* exp.(0.5*p) +0.5*p2 .* exp.(0.5*p) ) )*
(r[2]-r[1])
en = en+sum(0.25*exp.(2*p)-mu*exp.(p))*(r[2]-r[1])
en = en - 2.0*sum(exp.(p) .* g)*(r[2]-r[1]) en = en - 2.0*sum(exp.(p) .* g)*(r[2]-r[1])
nonloc = sum(exp.(p) .* (ps) .* exp.(p))*(r[2]-r[1])./(exp.(p)) nonloc = (sum(exp.(p) .* (ps) .* exp.(p))*(r[2]-r[1])) ./(exp.(p))
en = en - sum(nonloc)*(r[2]-r[1]) en = en + sum(nonloc)*(r[2]-r[1])
return en return en
end end
function gradient_ng_ln!(G,x)
p = psi(x)
ps = dpsi(x)
p2 = ddpsi(x)
for i=1:nfuns
bs = basisfunc(i)
dbs = dbasisfunc(i)
ddbs = ddbasisfunc(i)
#Gradient contribution arising from delta-interaction and mu
G[i] = sum(0.5*bs .* exp.(2*p) - mu * bs .* exp.(p))*(r[2]-r[1])
#Gradient contribution arising from the von Weizsäcker term
G[i] = G[i] - 0.5 * sum(0.5*bs .* ( 0.25*ps.^2.0 .* exp.(0.5*p) +0.5*p2 .* exp.(0.5*p)))*
(r[2]-r[1])
G[i] = G[i] - 0.5 * sum(exp.(0.5*p) .* (0.5*ps .* dbs .* exp.(0.5*p) +
0.25*0.5*ps.^2.0 .* bs .* exp.(0.5*p) + 0.5*ddbs .* exp.(0.5*p) +
0.25*p2 .* bs .* exp.(0.5*p)) )
#Gradient contribution arising from the nonlocal part
nonloc1 = (sum(2.0*bs .* exp.(2*p) .* ps )*(r[2]-r[1])) ./(exp.(p))
nonloc2 = (sum(exp.(2*p) .* dbs )*(r[2]-r[1])) ./(exp.(p))
nonloc3 = (sum(-exp.(p) .* (ps) .* exp.(p))*(r[2]-r[1])) ./(exp.(2*p))
G[i] = G[i] + sum(nonloc1 + nonloc2 + nonloc3 .* bs .*exp.(p) )
end
end
#td = TwiceDifferentiable(energy_tf,0.1*ones(60); autodiff=:forward) #td = TwiceDifferentiable(energy_tf,0.1*ones(60); autodiff=:forward)
res = optimize(energy_tf_vw,gradient_tf_vw!,0.1*ones(nfuns),LBFGS(),Optim.Options(iterations=10000)) initial = 0.1*zeros(nfuns)
initial[5] = 1.0
#initial[5] = 1.0
res = optimize(energy_ng,gradient_ng!,initial,LBFGS(),Optim.Options(iterations=50))
print(res) print(res)
finalrho = exp.(psi(Optim.minimizer(res))) finalrho = psi(Optim.minimizer(res)).^2.0
#finalrho = psi(Optim.minimizer(res)).^2.0 #finalrho = exp.(psi(Optim.minimizer(res)))
print(sum(finalrho)*(r[2]-r[1])) print(sum(finalrho)*(r[2]-r[1]))
plot(r,finalrho)
\ No newline at end of file plot(r,g)
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment