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

began slater impl

parent dec92943
Branches
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Makefile
+ 1
− 1
View file @ cc434609
OBJS = types.o constants.o data.o utils.o derivatives.o trapezoidalrule.o \
hartreepot.o mesh.o solverhovp.o thomasfermi.o \
hartreepot.o mesh.o solverhovp.o thomasfermi.o kedf_nagy.o minimize_slatergrid.o\
xc.o initialize.o energies.o energyminimizer.o main_Be.o
FC = gfortran
......
fit.log 0 → 100644
+ 51
− 0
View file @ cc434609
*******************************************************************************
Tue Jul 23 10:30:55 2019
FIT: data read from "results" using 1:3
format = x:z
#datapoints = 2001
residuals are weighted equally (unit weight)
function used for fitting: f(r)
f(r) = (sqrt(a**3/3.141592653) * exp(-a*r) * 4*3.141592653 * r**2)**2 + (sqrt(b**5/(3*3.141592653)) * r * exp(-b*r) * 4*3.141592653 * r**2)**2
fitted parameters initialized with current variable values
*******************************************************************************
Tue Jul 23 10:31:33 2019
FIT: data read from "results" using 1:3
format = x:z
#datapoints = 2001
residuals are weighted equally (unit weight)
function used for fitting: f(x)
f(x) = (sqrt(a**3/3.141592653) * exp(-a*x) * 4*3.141592653 * x**2)**2 + (sqrt(b**5/(3*3.141592653)) * x * exp(-b*x) * 4*3.141592653 * x**2)**2
fitted parameters initialized with current variable values
iter chisq delta/lim lambda a b
0 1.8987308053e+05 0.00e+00 1.05e+01 1.000000e+00 1.000000e+00
20 3.0950350096e+02 -1.08e-01 1.05e-13 1.106884e+01 8.776037e+00
After 20 iterations the fit converged.
final sum of squares of residuals : 309.504
rel. change during last iteration : -1.07622e-06
degrees of freedom (FIT_NDF) : 1999
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 0.393483
variance of residuals (reduced chisquare) = WSSR/ndf : 0.154829
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = 11.0688 +/- 0.1746 (1.577%)
b = 8.77604 +/- 0.05125 (0.584%)
correlation matrix of the fit parameters:
a b
a 1.000
b -0.169 1.000
kedf_nagy.f90 0 → 100644
+ 46
− 0
View file @ cc434609
module kedf_nagy
use types, only: dp
use derivatives, only: rDiff1, rDiff3
use data, only: Nmesh, r, Dr, Z, Norbs
use trapezoidalrule, only: trapz7
implicit none
public CalculateNagyPotential
contains
subroutine CalculateNagyPotential(rho, vks)
real(dp) :: rho(:), vks(:), dvksdr(Nmesh)
real(dp) :: energy, t(Nmesh), h(Nmesh), tpredict(Nmesh)
real(dp) :: k1,k2,k3,k4
integer :: i
end subroutine CalculateNagyPotential
subroutine tsolver(rho, vks, t)
real(dp) :: rho(:), vks(:)
real(dp) :: d3rhodr3(Nmesh), dvksdr(Nmesh), h(Nmesh)
real(dp) :: t(Nmesh), Tr
integer :: i
real(dp) :: k1, k2, k3, k4
t = 0.0_dp
d3rhodr3 = rDiff3(rho)
dvksdr = rDiff1(vks)
!We can solve kinetic energy density directly from Nagy's(2008) equation (22)
!Starting from initial conndition that kinetic energy density behavse correctly
!asymptotically.
h(:) = -Dr(:)
do i = 1680,2,-1
k1 = (h(i)) * ((-1.0_dp/8.0_dp) * d3rhodr3(i) - 0.5_dp * rho(i) * dvksdr(i))
k2 = (h(i) + h(i-1))*0.5_dp * ((-1.0_dp/8.0_dp) * (d3rhodr3(i) + d3rhodr3(i-1) )*0.5_dp -0.5_dp * (rho(i)*dvksdr(i) + &
rho(i-1)*dvksdr(i-1))*0.5_dp )
k3 = (h(i) + h(i-1))*0.5_dp * ((-1.0_dp/8.0_dp) * (d3rhodr3(i) + d3rhodr3(i-1) )*0.5_dp -0.5_dp * (rho(i)*dvksdr(i) + &
rho(i-1)*dvksdr(i-1))*0.5_dp )
k4 = (h(i)) * ((-1.0_dp/8.0_dp) * d3rhodr3(i-1) - 0.5_dp * rho(i-1) * dvksdr(i-1))
t(i-1) = t(i) + k1/6.0_dp + k4/6.0_dp + k2/3.0_dp + k3/3.0_dp
end do
!Tr = trapz7(t, 1, Nmesh)
!print *, Tr
end subroutine tsolver
end module kedf_nagy
\ No newline at end of file
kedf_nagy.mod 0 → 100644
+ 0
− 0
View file @ cc434609
File added
kedf_nagy.o 0 → 100644
+ 0
− 0
View file @ cc434609
File added
main_Be
+ 0
− 0
View file @ cc434609
No preview for this file type
main_Be.f90
+ 12
− 9
View file @ cc434609
......@@ -2,7 +2,7 @@ program main_Be
use types, only: dp
use constants, only: pi
use data
use energyminimizer, only: tsolver, ConjugateGradientMethod
use energyminimizer, only:ConjugateGradientMethod
use init, only: initialize
use derivatives, only: genDcoeffs,xDiff1,xDiff2,xDiff3, Weizsacker
use derivatives, only: rDiff1,rDiff2,rDiff3
......@@ -13,6 +13,8 @@ program main_Be
use solverhovp, only: ofdft_step
use utils, only: savetxt
use energies, only: get_energies
use kedf_nagy, only: CalculateNagyPotential, tsolver
use minimize_slatergrid, only: SlaterDensity
implicit none
integer :: i,tmp,iterMain
real(dp), allocatable :: f(:),dfdx(:),dfdr(:),d2fdr2(:),d3fdr3(:)
......@@ -76,9 +78,9 @@ program main_Be
vH=hartree(RhoOld)
vXC=LDAvX(RhoOld)+LDAvC(RhoOld)
vKS=vExt+vH+vXC
!
ConvTarget=1.0e-5
allocate(initrho(Nmesh))
! The main iteration loop can now begin
do iterMain=1,itersMainMax
print *,''
......@@ -111,6 +113,7 @@ program main_Be
kin_up = (1 + (4*3.141592653_dp)**2 * (Z/2.0_dp)**(2.0_dp/3.0_dp)*9.0_dp/5**(2.0_dp/3.0_dp) )*vw
rhomu = trapz7(muNew*RhoNew(:,1),1,Nmesh)
vptest = xclb - vw -gkinetic+rhomu
call tsolver(RhoNew(:,1),vks(:), initrho(:))
! Check convergence
if (L2Norm<ConvTarget) then
write(*,'(A48)') '################################################'
......@@ -126,6 +129,7 @@ program main_Be
write(*,'(A29,F11.6)') 'Electron-nucleus energy =',ENuc
write(*,'(A29,F11.6)') 'Exchange-correlation energy =',EXC
write(*,'(A40)') '----------------------------------------'
write(*,'(A29,F11.6)') 'test kinetic energy =',trapz7(initrho(:),1,Nmesh)
write(*,'(A29,F11.6)') 'Von W kinetic energy =',vw
write(*,'(A29, F11.6)')'Von W lower bound =',vw_lb
write(*,'(A29, F11.6)')'Von W upper bound =',vw_up
......@@ -154,6 +158,7 @@ program main_Be
write(*,'(A29,F11.6)') 'Electron-nucleus energy =',ENuc
write(*,'(A29,F11.6)') 'Exchange-correlation energy =',EXC
write(*,'(A40)') '----------------------------------------'
write(*,'(A29,F11.6)') 'test kinetic energy =',trapz7(initrho(:),1,Nmesh)
write(*,'(A29,F11.6)') 'Von W kinetic energy =',vw
write(*,'(A29, F11.6)')'Von W lower bound =',vw_lb
write(*,'(A29, F11.6)')'Von W upper bound =',vw_up
......@@ -179,11 +184,9 @@ program main_Be
allocate(vwp(Nmesh))
allocate(muar(Nmesh))
allocate(betaterm2(Nmesh))
allocate(initrho(Nmesh))
initrho = (8.0_dp*exp(-2.0_dp*R*Z)*R**2.0_dp*Z**3.0_dp+&
exp(-R*Z)*R**2.0_dp*Z**3.0_dp*(1.0_dp-(R*Z)/2.0_dp)**2.0_dp)
call ConjugateGradientMethod(RhoNew(:,1))
betaterm2(:) = rDiff1(RhoNew(:,1))
betaterm2 = SlaterDensity(5.0_dp, 1.3_dp)
print *, trapz7(betaterm2,1,Nmesh)
vwp = Weizsacker(RhoNew(:,1))
muar = muNew
outArray(:,1)=r
......@@ -195,7 +198,7 @@ program main_Be
!outArray(:,3)=vPNew(:,2)
outArray(:,3) = RhoNew(:,1)
!outArray(:,4)=vPNew(:,1)
outArray(:,2)=vPNew(:,1)
outArray(:,2)= betaterm2(:)
!outArray(:,4)= vpint(:)
!outArray(:,5)=
call savetxt(filename,outArray)
......
main_Be.o
+ 0
− 0
View file @ cc434609
No preview for this file type
minimize_slatergrid.f90 0 → 100644
+ 50
− 0
View file @ cc434609
module minimize_slatergrid
use types, only: dp
use trapezoidalrule, only: trapz7
use data, only: r, Dr, Z, Norbs, Nmesh
use xc, only: LDAeX, LDAvC
use derivatives, only: Weizsacker, rDiff1, rDiff2, rDiff3
use kedf_nagy, only: tsolver
implicit none
contains
recursive function factorial(n) result(f)
integer :: n, f
if(n == 1) then
f=1
else
f = n * factorial(n-1)
end if
end function factorial
function SlaterDensity(a,b) result(dens)
real(dp) :: dens(Nmesh)
real(dp) :: coeffs(Norbs)
real(dp) :: orb1(Nmesh), orb2(Nmesh), dens1(Nmesh), dens2(Nmesh)
real(dp) :: a,b
orb1(:) = sqrt(a**3/3.141592653) * exp(-a*r(:)) * 4*3.141592653 * r(:)**2
orb2(:) = sqrt(b**5/(3*3.141592653)) * r(:) * exp(-b*r(:)) * 4*3.141592653 * r(:)**2
dens1(:) = 2/trapz7(orb1(:)**2, 1, Nmesh) * orb1(:)**2
dens2(:) = 2/trapz7(orb2(:)**2, 1, Nmesh) * orb2(:)**2
dens(:) = dens1(:) + dens2(:)
end function SlaterDensity
subroutine GenerateSlaterGrid(npoints, acoeffs, bcoeffs)
integer :: npoints, i
real(dp) :: acoeffs(:), bcoeffs(:)
do i=1, npoints
acoeffs(i) = 3.5 + i*0.05
bcoeffs(i) = 0.3 + i*0.05
end do
end subroutine GenerateSlaterGrid
end module minimize_slatergrid
\ No newline at end of file
minimize_slatergrid.mod 0 → 100644
+ 0
− 0
View file @ cc434609
File added
minimize_slatergrid.o 0 → 100644
+ 0
− 0
View file @ cc434609
File added
results
+ 2001
− 2001
View file @ cc434609
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment