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

attempted implementation of another functional

parent d8391535
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
basis.f90
+ 1
− 1
View file @ 9cc88f79
...@@ -32,7 +32,7 @@ function create_basis(nelectrons,nfuns,ngrid,L) result(bas) ...@@ -32,7 +32,7 @@ function create_basis(nelectrons,nfuns,ngrid,L) result(bas)
x(i) = (i) * ep x(i) = (i) * ep
end do end do
do i=1, nfuns do i=1, nfuns
f(:,i) = sqrt(2/L) * sin(i*pi*x(:)/L) f(:,i) = sqrt(2/L) * sin(i*pi*x(:)/L)!+1.0_dp
fder(:,i) = sqrt(2/L) * i*pi/L *cos(i*pi*x(:)/L) fder(:,i) = sqrt(2/L) * i*pi/L *cos(i*pi*x(:)/L)
end do end do
bas%nfuns = nfuns bas%nfuns = nfuns
......
energies.f90
+ 23
− 59
View file @ 9cc88f79
...@@ -55,28 +55,7 @@ function get_exchange_energy(bas,nx) result(ex) ...@@ -55,28 +55,7 @@ function get_exchange_energy(bas,nx) result(ex)
end function end function
function get_tf_pot(bas,alpha,sn,vec,mu) result(k)
type(basis_t) :: bas
real(dp) :: k,sn(:),vec(:),alpha,beta, kindens(size(bas%x)), deldens(size(bas%x)),mu,const
real(dp) :: psi(size(bas%x)), pi
pi = 4.0_dp*atan(1.0_dp)
psi = expand_psi_ofdft(bas,vec(:)+alpha*sn(:))
k = sum( abs(2.0_dp*psi(:)*(0.5_dp * psi(:)**2.0_dp+ pi**2.0_dp/2.0_dp *psi(:)**4.0_dp-mu) )**2.0_dp*bas%dx)
!print *, k
end function
function get_tf_pot_grad(bas,vec,mu) result(grad)
type(basis_t) :: bas
real(dp) :: grad(bas%nfuns), vec(:), mu, psi(size(bas%x)),pi
integer :: i
pi = 4.0_dp*atan(1.0_dp)
psi = expand_psi_ofdft(bas,vec)
grad = 0.0_dp
do i=1, bas%nfuns
grad(i) = grad(i)+sum(bas%dx*(4.0_dp*bas%f(:,i)*(0.5_dp * psi(:)**2.0_dp+ pi**2.0_dp/2.0_dp *psi(:)**4.0_dp-mu)))
grad(i) = grad(i)+sum(bas%dx*(4.0_dp*psi(:)*(bas%f(:,i)*psi(:)+2.0_dp*pi**2.0_dp *bas%f(:,i)*psi(:)**3.0_dp)) )
end do
end function
function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k) function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k)
type(basis_t) :: bas type(basis_t) :: bas
real(dp) :: k, g(size(bas%x)), beta(size(bas%x)),pi,alpha,mu,rho(size(bas%x)), drho(size(bas%x)) real(dp) :: k, g(size(bas%x)), beta(size(bas%x)),pi,alpha,mu,rho(size(bas%x)), drho(size(bas%x))
...@@ -104,38 +83,31 @@ function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k) ...@@ -104,38 +83,31 @@ function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k)
g(:) = beta(:)/rho(:) g(:) = beta(:)/rho(:)
drho(:) = derivative(bas,rho) drho(:) = derivative(bas,rho)
nonloc(:) = 2*g(:)*drho(:) nonloc(:) = sum(derivative(bas,psi(:)**2.0_dp)*psi(:)**2.0_dp * bas%dx)/psi(:)**2.0_dp
!goto 10 !goto 10
do i=size(bas%x)/2+1, size(bas%x) goto 20
do j=1, size(bas%x)-i+1
xint(i) = xint(i) - g(i+j-1)*drho(i+j-1)*psi(i+j-1)**4.0_dp*dx
!xint(i) = xint(i) - 2.0_dp*g(i+j-1)*drho(i+j-1)*rho(i+j-1)**2.0_dp
end do
end do
xint(size(bas%x)/2) = xint(size(bas%x)/2+1)
do i=1, size(bas%x)/2-1
xint(i) = xint( size(bas%x)-i+1)
end do
!10 continue
xint(:) = xint(:)/(rho**2.0_dp)
!goto 20
open(12, file="gnonloc.dt") open(12, file="gnonloc.dt")
do i=1,size(bas%x) do i=1,size(bas%x)
write(12,*) bas%x(i), xint(i),g(i),drho(i),psi(i)**4 write(12,*) bas%x(i), nonloc(i),g(i),drho(i),psi(i)**2.0_dp
end do end do
close(12) close(12)
!20 continue stop
20 continue
k = sum(( abs( derivative(bas,psi(:)**2.0_dp) ) )**2.0_dp/(psi(:)**2.0_dp)*dx)*1.0_dp/(8.0_dp) k = sum(( abs( derivative(bas,psi(:)**2.0_dp) ) )**2.0_dp/(psi(:)**2.0_dp)*dx)*1.0_dp/(8.0_dp)
k = k+mupar*sum(psi(:)**2.0_dp*dx) - 2.0_dp * sum(g(:)*psi(:)**2.0_dp*dx) & k = k+mupar*sum(psi(:)**2.0_dp*dx) - 2.0_dp * sum(g(:)*psi(:)**2.0_dp*dx) &
+sum(0.25_dp*psi(:)**4.0_dp*dx)-mu*sum(psi(:)**2.0_dp*dx)+ sum(xint(:)*dx) +sum(0.25_dp*psi(:)**4.0_dp*dx)-mu*sum(psi(:)**2.0_dp*dx)-sum(nonloc(:)*bas%dx)
!print *, "here"
!print *, k !print *, k
!stop
end function end function
function get_nagy_gradient(bas,vec,mu) result(grad) function get_nagy_gradient(bas,vec,mu) result(grad)
type(basis_t) :: bas type(basis_t) :: bas
real(dp) :: vec(:),mu,psi(size(bas%x)),grad(bas%nfuns),xint(size(bas%x)),drho(size(bas%x)),mupar real(dp) :: vec(:),mu,psi(size(bas%x)),grad(bas%nfuns),xint(size(bas%x)),drho(size(bas%x)),mupar
real(dp) :: nonloc(size(bas%x)),rho(size(bas%x)),beta(size(bas%x)),dx,g(size(bas%x)),pi real(dp) :: nonloc(size(bas%x),bas%nfuns),rho(size(bas%x)),beta(size(bas%x)),dx,g(size(bas%x)),pi
integer :: i,j,k,z integer :: i,j,k,z
!Beta term for non-interacting system !Beta term for non-interacting system
beta = 0.0_dp beta = 0.0_dp
...@@ -160,26 +132,17 @@ function get_nagy_gradient(bas,vec,mu) result(grad) ...@@ -160,26 +132,17 @@ function get_nagy_gradient(bas,vec,mu) result(grad)
g(:) = beta(:)/rho(:) g(:) = beta(:)/rho(:)
drho(:) = derivative(bas,rho) drho(:) = derivative(bas,rho)
nonloc(:) = 2*g(:)*drho(:) nonloc = 0.0_dp
do i=1, bas%nfuns
nonloc(:,i) = nonloc(:,i)+ 2.0_dp*sum((bas%fder(:,i)*psi(:)+&
bas%f(:,i)*derivative(bas,psi(:)) )*psi(:)**2.0_dp*bas%dx)/psi(:)**2.0_dp
nonloc(:,i) = nonloc(:,i) + sum(bas%f(:,i)*derivative(bas,psi(:)**2.0_dp)*bas%dx)/psi(:)**2.0_dp
nonloc(:,i) = nonloc(:,i) - 2.0_dp*bas%f(:,i)*sum(psi(:)**2.0_dp*derivative(bas,psi(:)**2.0_dp)*bas%dx)/psi(:)**3.0_dp
end do
!grad(:) = get_vw_gradient !grad(:) = get_vw_gradient
!goto 10 !goto 10
do k=1, bas%nfuns
do i=size(bas%x)/2+1, size(bas%x)
do j=1, size(bas%x)-i+1
!print *, i
xint(i) = xint(i) - 4.0_dp*g(i+j-1)*drho(i+j-1)*bas%f(i+j-1,k)*psi(i+j-1)**3.0_dp*dx
end do
end do
xint(size(bas%x)/2) = xint(size(bas%x)/2+1)
do z=1,size(bas%x)/2-1
xint(z) = xint( size(bas%x)-z+1)
end do
grad(k) = grad(k) + sum(xint(:)/(rho(:)**2.0_dp) *dx)
xint = 0.0_dp
!print *, k !10 continu
end do
!10 continue
do i=1, bas%nfuns do i=1, bas%nfuns
!print *, i !print *, i
grad(i) = grad(i) + sum(2.0_dp*mupar * bas%f(:,i) *psi(:)* dx)-sum(4.0_dp*g(:)*bas%f(:,i)*psi(:)*dx) grad(i) = grad(i) + sum(2.0_dp*mupar * bas%f(:,i) *psi(:)* dx)-sum(4.0_dp*g(:)*bas%f(:,i)*psi(:)*dx)
...@@ -191,6 +154,7 @@ function get_nagy_gradient(bas,vec,mu) result(grad) ...@@ -191,6 +154,7 @@ function get_nagy_gradient(bas,vec,mu) result(grad)
grad(i) = grad(i) - 2.0_dp/(8.0_dp) *& grad(i) = grad(i) - 2.0_dp/(8.0_dp) *&
sum(bas%dx * bas%f(:,i)*abs(derivative(bas,psi**2.0_dp))**2.0_dp/(psi**3.0_dp)) sum(bas%dx * bas%f(:,i)*abs(derivative(bas,psi**2.0_dp))**2.0_dp/(psi**3.0_dp))
grad(i) = grad(i) - sum(nonloc(:,i)*bas%dx)
end do end do
!print *, "we here?" !print *, "we here?"
end function end function
......
minim
+ 0
− 0
View file @ 9cc88f79
No preview for this file type
minim.f90
+ 7
− 6
View file @ 9cc88f79
...@@ -33,12 +33,13 @@ L = 8.0_dp ...@@ -33,12 +33,13 @@ L = 8.0_dp
!g_ptr => get_vw_gradient !g_ptr => get_vw_gradient
!f_ptr =>get_nagy_total_energy !f_ptr =>get_nagy_total_energy
!g_ptr => get_nagy_gradient !g_ptr => get_nagy_gradient
!f_ptr => get_vw_total_energy f_ptr => get_vw_total_energy
!g_ptr => get_vw_gradient g_ptr => get_vw_gradient
f_ptr => get_tf_pot !f_ptr => get_tf_pot
g_ptr => get_tf_pot_grad !g_ptr => get_tf_pot_grad
bas = basis_t(nel,20,2000,L) bas = basis_t(nel,60,2000,L)
!f_ptr => get_thomas_fermi_total_energy
!g_ptr => get_thomas_fermi_gradient
!call deq_optimize(bas,f_ptr,g_ptr) !call deq_optimize(bas,f_ptr,g_ptr)
call optimizer(bas,f_ptr,g_ptr) call optimizer(bas,f_ptr,g_ptr)
!L = get_nagy_total_energy(bas,0.0_dp,0.0_dp) !L = get_nagy_total_energy(bas,0.0_dp,0.0_dp)
......
optimize.f90
+ 14
− 7
View file @ 9cc88f79
...@@ -56,10 +56,11 @@ subroutine optimizer(bas,defunc,gradfunc) ...@@ -56,10 +56,11 @@ subroutine optimizer(bas,defunc,gradfunc)
procedure(func) :: defunc procedure(func) :: defunc
procedure(gs) :: gradfunc procedure(gs) :: gradfunc
real(dp) :: mu, a,b,c, optres real(dp) :: mu, a,b,c, optres
a = -50.0_dp a = 0.0_dp
b = 50.0_dp b = 50.0_dp
c = 25.0_dp c = 25.0_dp
optres = 10.0_dp optres = 10.0_dp
!goto 10
do while(abs(optres) > 1e-2_dp) do while(abs(optres) > 1e-2_dp)
c = (a+b)/2.0_dp c = (a+b)/2.0_dp
optres = deq_optimize(bas,defunc,gradfunc,c) optres = deq_optimize(bas,defunc,gradfunc,c)
...@@ -71,23 +72,25 @@ subroutine optimizer(bas,defunc,gradfunc) ...@@ -71,23 +72,25 @@ subroutine optimizer(bas,defunc,gradfunc)
print *, a,b print *, a,b
end do end do
print *,"mu now", c print *,"mu now", c
!10 continue
!optres = deq_optimize(bas,defunc,gradfunc,a)
end subroutine end subroutine
function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
type(Basis_t) :: bas type(Basis_t) :: bas
real(dp) :: grad(bas%nfuns),gradprev(bas%nfuns),b,stp,etp real(dp) :: grad(bas%nfuns),gradprev(bas%nfuns),b,stp,etp
real(dp) :: sn(bas%nfuns),snprev(bas%nfuns), curr(bas%nfuns),prev(bas%nfuns),gamma,mu, psi(size(bas%x)) real(dp) :: sn(bas%nfuns),snprev(bas%nfuns), curr(bas%nfuns),prev(bas%nfuns),gamma,mu, psi(size(bas%x))
real(dp) :: rho(size(bas%x)), prevenergy,currenergy,nelectron real(dp) :: rho(size(bas%x)), prevenergy,currenergy,nelectron,pi
integer :: i,j integer :: i,j
procedure(func) :: defunc procedure(func) :: defunc
procedure(gs) :: gradfunc procedure(gs) :: gradfunc
!mu = 8.7_dp !mu = 8.7_dp
pi = 4.0_dp*atan(1.0_dp)
stp = 0.0_dp stp = 0.0_dp
etp = 10.0_dp etp = 10.0_dp
prev(:) = 0.0_dp prev(:) = 0.0_dp
prev(:) = 0.01_dp prev(5) = 0.1_dp
!pi = defunc(bas,0.0_dp,prev,prev,mu)
psi(:) = expand_psi_ofdft(bas,prev) psi(:) = expand_psi_ofdft(bas,prev)
!print *, (bas%x(2)-bas%x(1))*sum(psi(:)**2.0_dp) !print *, (bas%x(2)-bas%x(1))*sum(psi(:)**2.0_dp)
curr(:) = 0.0_dp curr(:) = 0.0_dp
...@@ -104,6 +107,9 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) ...@@ -104,6 +107,9 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
prevenergy = defunc(bas,stp,sn(:),curr(:),mu) prevenergy = defunc(bas,stp,sn(:),curr(:),mu)
do i=1, 100000 do i=1, 100000
psi(:) = expand_psi_ofdft(bas,curr)
!print *, sum(psi(:)),"lol"
!print *, "mi"
gradprev(:) = grad(:) gradprev(:) = grad(:)
grad(:) = gradfunc(bas,curr,mu) grad(:) = gradfunc(bas,curr,mu)
b = beta(grad,gradprev) b = beta(grad,gradprev)
...@@ -126,7 +132,7 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) ...@@ -126,7 +132,7 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
print *, i, sum(grad*grad) print *, i, sum(grad*grad)
if(sum(grad(:)*grad(:)) < 1e-8_dp .or. & if(sum(grad(:)*grad(:)) < 1e-8_dp .or. &
(abs(defunc(bas,stp,sn(:),curr(:),mu)-prevenergy )<1e-8_dp .and.sum(grad(:)*grad(:)) < 2.0_dp ) ) then (abs(defunc(bas,stp,sn(:),curr(:),mu)-prevenergy )<1e-8_dp .and.sum(grad(:)*grad(:)) < 2.0_dp ) ) then
print *, "stopped",i !print *, "stopped",i
...@@ -134,10 +140,11 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) ...@@ -134,10 +140,11 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
psi = expand_psi_ofdft(bas,curr(:)) psi = expand_psi_ofdft(bas,curr(:))
rho = psi**2.0_dp rho = psi**2.0_dp
nelectron= sum(rho)*(bas%x(2)-bas%x(1))-bas%nelectrons nelectron= sum(rho)*(bas%x(2)-bas%x(1))-bas%nelectrons
print *, nelectron, defunc(bas,0.0_dp,sn,curr,mu)+mu*sum(bas%dx*psi**2) print *, "electron difference",nelectron
do j=1, size(bas%x) do j=1, size(bas%x)
write(12, *) bas%x(j), rho(j) write(12, *) bas%x(j), rho(j)
end do end do
!print *, curr
close(12) close(12)
exit exit
end if end if
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment