diff --git a/analyze_vp.py b/analyze_vp.py
new file mode 100644
index 0000000000000000000000000000000000000000..2b3126c0f5c67137b11c94a922420bfcc24070b0
--- /dev/null
+++ b/analyze_vp.py
@@ -0,0 +1,17 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+vpcont = np.loadtxt("cont-vp.dt")
+vpr = np.loadtxt("vp4.dt")
+
+x = np.linspace(0,8,2000)
+plt.rcParams.update({'font.size': 15})
+plt.xlabel("r (Bohr)")
+plt.ylabel("Pauli potential (a.u)")
+plt.plot(x,vpr[:,1],label="Correct Pauli potential")
+plt.plot(x,vpcont[:,1],label="Approximated Pauli potential")
+#plt.plot(x,rhos[:,3],label="Corrected potential")
+#plt.plot(betar[:,0],betar[:,1],label="Tarkka")
+plt.legend()
+plt.show()
\ No newline at end of file
diff --git a/analyze_vp_error.py b/analyze_vp_error.py
index 4df0dfbb416dad428600687821cb9a3abc3f8161..e54afc8e3b0bd062e82742d21fa0ced0ad1148eb 100644
--- a/analyze_vp_error.py
+++ b/analyze_vp_error.py
@@ -10,6 +10,7 @@ plt.xlabel("r (Bohr)")
plt.ylabel("Density (a.u)")
plt.plot(x,rhos[:,1],label="Correct density")
plt.plot(x,rhos[:,2],label="Slighty wrong potential")
+plt.plot(x,rhos[:,3],label="Corrected potential")
#plt.plot(betar[:,0],betar[:,1],label="Tarkka")
plt.legend()
plt.show()
\ No newline at end of file
diff --git a/hamilton.f90 b/hamilton.f90
index 953b5c46a176d71c8355318a8201f049c9b5eedc..5137d6f0f2f1c312cb3e60dfd20a54cf274cc1e4 100644
--- a/hamilton.f90
+++ b/hamilton.f90
@@ -2,7 +2,7 @@ module hamilton
use types, only: dp
use basis, only: basis_t, derivative
implicit none
-
+public
contains
function fnm(x,L,n,m) result(res)
@@ -498,33 +498,70 @@ subroutine get_pauli_ready(H,bas,vp)
end do
end subroutine
-subroutine get_pauli_continuous_appro(H,bas,nx)
+subroutine get_pauli_continuous_appro(H,bas,nx,g)
real(dp), intent(inout) :: H(:,:)
type(basis_t) :: bas
- real(dp) :: nx(:),vp(size(nx)),f(size(nx)), g(size(nx)), nxsq(size(nx))
- real(dp) :: dnxsq(size(nx)),ddnxsq(size(nx)), dnx(size(nx))
- real(dp) :: dx
+ real(dp) :: nx(:),vp(size(nx)),f(size(nx)), g(:), nxsq(size(nx))
+ real(dp) :: dnxsq(size(nx)),ddnxsq(size(nx)), dnx(size(nx)),pi,mu,c,nonloc(size(nx))
integer :: i,j
- dx = bas%x(2)-bas%x(1)
- nxsq(:) = sqrt(nx(:))
- dnxsq(:) = derivative(bas,nxsq(:))
- ddnxsq(:) = derivative(bas,dnxsq(:))
- dnx(:) = derivative(bas,nx(:))
- g(:) = -1/sqrt(nx) * 1.0_dp/2.0_dp *ddnxsq(:) - nx(:)**2.0_dp * (4.0_dp*atan(1.0_dp) )**2.0_dp/6.0_dp + 8.97_dp
- g(:) = 0.5_dp *g(:)
- f(:) = g(:) * derivative(bas,nx)**2.0_dp
+ pi = 4.0_dp * atan(1.0_dp)
+ c = 0.5_dp
+ mu = 5.0_dp**2.0_dp*pi**2.0_dp/(2.0_dp*8.0_dp**2.0_dp)
+ dnx = derivative(bas,nx)
+ f(:) =c *g(:)* dnx(:) *nx(:)
vp(:) = 0.0_dp
do i=1, size(bas%x)
do j=1, size(bas%x)-i+1
- vp(i) = vp(i)- f(i+j-1)*dx
+ vp(i) = vp(i)- f(i+j-1)*bas%dx
end do
end do
- vp(:) = -2.0_dp/nx(:)**2.0_dp * vp(:) !+ 2*g(:)
+ nonloc(:) = 2.0_dp/nx(:)**2.0_dp * vp(:)
+ vp(:) = 2.0_dp/nx(:)**2.0_dp * vp(:) - (2.0_dp-c)*g(:) +mu
+ print *, mu
open(12, file="cont-vp.dt")
do i=1, size(nx(:))
- write(12, * ) bas%x(i), vp(i), 2*g(i)*dnx(i)/nx(i)**2
+ write(12, * ) bas%x(i), vp(i), g(i), nonloc(i)
end do
close(12)
+end subroutine
+subroutine get_pauli_exact2_symmetry(H,bas,mu,beta,nx)
+ real(dp),intent(inout) :: H(:,:)
+ real(dp),intent(in) :: mu,nx(:)
+ type(basis_t),intent(in):: bas
+ real(dp) :: beta(:), vp(size(nx)), f(size(nx)),dnx(size(nx))
+ integer :: i,j,k
+ vp(:) = mu
+ dnx(:) = derivative(bas,nx(:))
+
+ vp(:) = vp(:) -2.0_dp*beta(:)/(nx(:))
+ f(:) = 0.0_dp
+ do i=1, size(bas%x)
+ do j=1, size(bas%x)-i+1
+ f(i) = f(i)- dnx(i+j-1)*beta(i+j-1)*bas%dx
+ end do
+ end do
+
+
+ vp(:) = vp(:) + 2.0_dp*f(:)/(nx(:))**2.0_dp
+
+
+ do i=1,999
+ vp(i+1001) = vp(1000-i)
+ end do
+ vp(1001) = vp(1000)
+ open(12, file="brho.dt")
+ do i=1, size(nx(:))
+ write(12, * ) bas%x(i), vp(i)
+ end do
+ close(12)
+ do i=1, size(bas%f(1,:))
+ do j=i, size(bas%f(1,:))
+ H(i,j) = H(i,j)+ get_matrix_element(bas,vp(:), i,j)
+ H(j,i) = H(i,j)
+ end do
+ end do
+
+
end subroutine
subroutine get_thomas_fermi(H,bas,nx)
type(basis_t) :: bas
diff --git a/main.f90 b/main.f90
index bf70ef427b823f2a69d599baa0f348080e34045a..ce2b390e8af98893b278bdd43033d03c9a1f6689 100644
--- a/main.f90
+++ b/main.f90
@@ -9,7 +9,7 @@ integer :: npeak
type(eigres_t) :: sol(7), sol2
type(basis_t) :: testb
type(dft_t) :: res, res1(5)
-real(dp):: H(1,100,100),a,b,L,V0,nx(2000),hart(5,2000), pi, psi(2000,5),nx3(2000),db,nx2(2000)
+real(dp):: H(1,100,100),a,b,L,V0,nx(2000),hart(2000), pi, psi(2000,5),nx3(2000),db,nx2(2000), g(2000)
real(dp):: ergs(5)
integer :: i,j,k,norbs,nelectrons
norbs = 5
@@ -25,36 +25,41 @@ V0 = 100.0_dp
H(:,:,:) = 0.0_dp
testb = basis_t(nelectrons,100,2000,L)
nx(:) = 0.0_dp
-
+hart(:) = 0.0_dp
do i=1, norbs
nx(:) = nx(:)+2.0_dp * 2.0_dp/testb%L * sin(i*pi/testb%L * testb%x(:))**2.0_dp
psi(:,i) = sqrt(2.0_dp/testb%L) * sin(i*pi/testb%L * testb%x(:))
ergs(i) = i**2.0_dp * pi**2.0_dp/(2.0_dp * L**2.0_dp)
+ hart(:) = hart(:) + 2.0_dp * ergs(i) * psi(:,i)**2.0_dp
end do
+g(:) = hart(:)/nx(:)
-!res1(2) = solve_dft(testb,nx,hart(1,:),.false.,norbs,.true.,10)
+call get_pauli_exact2_symmetry(H(1,:,:),testb,ergs(5),hart,nx)
+res1(2) = solve_dft(testb,nx,hart(:),.false.,norbs,.true.,10)
!print *, sum((res1(1)%nx(:) -res1(2)%nx(:))**2*(testb%x(2)-testb%x(1)) )
!print *, res1(2)%energies
!call get_pauli_beta_rho(H,testb,sol2%energies(2),db,rhobeta)
!call get_pauli_exact(H(1,:,:),testb,res1(2)%psi,res1(2)%energies,res1(2)%nx)
-call get_pauli_exact(H(1,:,:),testb,psi,ergs,nx)
-call get_exchange(H(1,:,:),testb,nx)
-call get_hartree(H(1,:,:),testb,nx)
-call get_zero_energy(H(1,:,:),testb)
-sol2 = eigres_t(1,size(testb%f(1,:)))
-sol2 = solve_eigenproblem(H(1,:,:),testb,1)
-print *, "hey"
-nx2 = expand_rho(testb,sol2%vec,.true.,1)
+!call get_pauli_exact(H(1,:,:),testb,psi,ergs,nx)
+!call get_exchange(H(1,:,:),testb,nx)
+!call get_hartree(H(1,:,:),testb,nx)
+!call get_zero_energy(H(1,:,:),testb)
+!sol2 = eigres_t(1,size(testb%f(1,:)))
+!sol2 = solve_eigenproblem(H(1,:,:),testb,1)
+!print *, "hey"
+!nx2 = expand_rho(testb,sol2%vec,.true.,1)
H = 0.0_dp
-call get_pauli_exact_symmetry(H(1,:,:),testb,psi,ergs,nx)
-call get_exchange(H(1,:,:),testb,nx)
-call get_hartree(H(1,:,:),testb,nx)
-call get_zero_energy(H(1,:,:),testb)
-sol2 = eigres_t(1,size(testb%f(1,:)))
-sol2 = solve_eigenproblem(H(1,:,:),testb,1)
-nx3 = expand_rho(testb,sol2%vec,.true.,1)
+call get_pauli_exact_symmetry(H(1,:,:),testb,res1(2)%psi,res1(2)%energies,res1(2)%nx)
+!
+!call get_exchange(H(1,:,:),testb,nx)
+!call get_hartree(H(1,:,:),testb,nx)
+!call get_zero_energy(H(1,:,:),testb)
+!sol2 = eigres_t(1,size(testb%f(1,:)))
+!sol2 = solve_eigenproblem(H(1,:,:),testb,1)
+!nx3 = expand_rho(testb,sol2%vec,.true.,1)
+call get_pauli_continuous_appro(H(1,:,:),testb,res1(2)%nx,g)
!call get_pauli_continuous_appro(H(1,:,:),testb,res1(2)%nx)
!print*, "Kinetic", get_kinetic_energy(testb,res1(2)%psi)
!print*, "Exchange", get_exchange_energy(testb,res1(2)%nx)
@@ -64,7 +69,7 @@ nx3 = expand_rho(testb,sol2%vec,.true.,1)
!print *, res1(2)%energies
open(12, file="rho.dt")
do i=1, size(testb%x)
- write(12, *) testb%x(i),nx(i),nx2(i),nx3(i)
+ write(12, *) testb%x(i),nx(i),res1(2)%nx(i)
end do
close(12)
diff --git a/solver.f90 b/solver.f90
index 67784cfa031fc9c67ab6ebae437e7a7174068246..032dfa17b895b27f1aa7e092a7a5f3e102ab57b9 100644
--- a/solver.f90
+++ b/solver.f90
@@ -139,7 +139,7 @@ function solve_dft(bas,initrho,beta,of,norbs,diisflag,ndiis) result(res)
dx = bas%x(2)-bas%x(1)
preve = -10000.0_dp
H = 0.0_dp
- alpha = 0.01_dp
+ alpha = 0.1_dp
ediff = 0.0_dp
sol = eigres_t(norbs,size(bas%f(1,:)))
sol%vec(:,:) = 0.0_dp