diff --git a/basis.f90 b/basis.f90 index de1400e4937e885ed082c7c514f1fe4af3bc59e2..bc4a8ae8bab844290754158d4489101cc151b0df 100644 --- a/basis.f90 +++ b/basis.f90 @@ -32,7 +32,7 @@ function create_basis(nelectrons,nfuns,ngrid,L) result(bas) x(i) = (i) * ep end do 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) end do bas%nfuns = nfuns diff --git a/energies.f90 b/energies.f90 index 8571d0a877c04a0d5c811a4594739f74b3ecbf09..a8a1a81d6a3ed075fcb8eedbd98318cbd8420c9e 100644 --- a/energies.f90 +++ b/energies.f90 @@ -55,28 +55,7 @@ function get_exchange_energy(bas,nx) result(ex) 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) 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)) @@ -104,38 +83,31 @@ function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k) g(:) = beta(:)/rho(:) drho(:) = derivative(bas,rho) - nonloc(:) = 2*g(:)*drho(:) - !goto 10 + nonloc(:) = sum(derivative(bas,psi(:)**2.0_dp)*psi(:)**2.0_dp * bas%dx)/psi(:)**2.0_dp - do i=size(bas%x)/2+1, size(bas%x) - 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 + + !goto 10 + + goto 20 open(12, file="gnonloc.dt") 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 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 = 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) - !print *, k + +sum(0.25_dp*psi(:)**4.0_dp*dx)-mu*sum(psi(:)**2.0_dp*dx)-sum(nonloc(:)*bas%dx) + !print *, "here" + !print *, k + !stop end function function get_nagy_gradient(bas,vec,mu) result(grad) 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) :: 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 !Beta term for non-interacting system beta = 0.0_dp @@ -160,26 +132,17 @@ function get_nagy_gradient(bas,vec,mu) result(grad) g(:) = beta(:)/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 !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 - end do - !10 continue + + !10 continu do i=1, bas%nfuns !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) @@ -191,6 +154,7 @@ function get_nagy_gradient(bas,vec,mu) result(grad) 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)) + grad(i) = grad(i) - sum(nonloc(:,i)*bas%dx) end do !print *, "we here?" end function diff --git a/minim b/minim index 5b3f8192f8a88289efbd7abc5a889dc21a92bf31..fd105a2c92a3c97ec465f96acd9ed19a9f6e8084 100755 Binary files a/minim and b/minim differ diff --git a/minim.f90 b/minim.f90 index a3edb238b6d24734ee640b85c045eae78823e263..538d4aa86ee7f9bb314b633ce68e45599ed8ec81 100644 --- a/minim.f90 +++ b/minim.f90 @@ -33,12 +33,13 @@ L = 8.0_dp !g_ptr => get_vw_gradient !f_ptr =>get_nagy_total_energy !g_ptr => get_nagy_gradient -!f_ptr => get_vw_total_energy -!g_ptr => get_vw_gradient -f_ptr => get_tf_pot -g_ptr => get_tf_pot_grad -bas = basis_t(nel,20,2000,L) - +f_ptr => get_vw_total_energy +g_ptr => get_vw_gradient +!f_ptr => get_tf_pot +!g_ptr => get_tf_pot_grad +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 optimizer(bas,f_ptr,g_ptr) !L = get_nagy_total_energy(bas,0.0_dp,0.0_dp) diff --git a/optimize.f90 b/optimize.f90 index d33b8ae2246327ec7b509dd4307e3cc9318aae9f..72db14efc16a421568dd8aa65359d55418f22166 100644 --- a/optimize.f90 +++ b/optimize.f90 @@ -56,10 +56,11 @@ subroutine optimizer(bas,defunc,gradfunc) procedure(func) :: defunc procedure(gs) :: gradfunc real(dp) :: mu, a,b,c, optres - a = -50.0_dp + a = 0.0_dp b = 50.0_dp c = 25.0_dp optres = 10.0_dp + !goto 10 do while(abs(optres) > 1e-2_dp) c = (a+b)/2.0_dp optres = deq_optimize(bas,defunc,gradfunc,c) @@ -71,23 +72,25 @@ subroutine optimizer(bas,defunc,gradfunc) print *, a,b end do print *,"mu now", c - + !10 continue + !optres = deq_optimize(bas,defunc,gradfunc,a) end subroutine function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) type(Basis_t) :: bas 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) :: rho(size(bas%x)), prevenergy,currenergy,nelectron + real(dp) :: rho(size(bas%x)), prevenergy,currenergy,nelectron,pi integer :: i,j procedure(func) :: defunc procedure(gs) :: gradfunc !mu = 8.7_dp + pi = 4.0_dp*atan(1.0_dp) stp = 0.0_dp etp = 10.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) !print *, (bas%x(2)-bas%x(1))*sum(psi(:)**2.0_dp) curr(:) = 0.0_dp @@ -104,6 +107,9 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) prevenergy = defunc(bas,stp,sn(:),curr(:),mu) do i=1, 100000 + psi(:) = expand_psi_ofdft(bas,curr) + !print *, sum(psi(:)),"lol" + !print *, "mi" gradprev(:) = grad(:) grad(:) = gradfunc(bas,curr,mu) b = beta(grad,gradprev) @@ -126,7 +132,7 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) print *, i, sum(grad*grad) 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 - print *, "stopped",i + !print *, "stopped",i @@ -134,10 +140,11 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron) psi = expand_psi_ofdft(bas,curr(:)) rho = psi**2.0_dp 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) write(12, *) bas%x(j), rho(j) end do + !print *, curr close(12) exit end if