continued work on paper; optimizer

analyze_vp_error.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+
+rhos = np.loadtxt("rho.dt")
+
+x = np.linspace(0,8,2000)
+plt.rcParams.update({'font.size': 15})
+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(betar[:,0],betar[:,1],label="Tarkka")
+plt.legend()
+plt.show()
\ No newline at end of file

hamilton.f90

@@ -297,7 +297,60 @@ subroutine get_pauli_exact(H,bas,psi,energies,nx)
     close(12)
     
 end subroutine
+subroutine get_pauli_exact_symmetry(H,bas,psi,energies,nx)
+    type(basis_t),intent(in):: bas
+    real(dp),intent(inout)  :: H(:,:)
+    real(dp),intent(in)     :: nx(:),energies(:),psi(:,:)
+    real(dp)                :: vp(size(bas%x)),vp2(size(bas%x)), f(size(bas%x))
+    real(dp)                :: pi,dx, eqenergies(bas%nelectrons/2),beta(size(bas%x))
+    real(dp)                :: dbeta(size(bas%x)), dnx(size(bas%x))
+    integer                 :: i,j
+    pi = 4.0_dp*atan(1.0_dp)
+    dx = bas%x(2)-bas%x(1)
+    vp(:) = 0.0_dp
+    f(:) = 0.0_dp
+    beta(:) = 0.0_dp
+    dnx = derivative(bas,nx(:))
+
+    do i=1, bas%nelectrons/2
+        beta(:) = beta(:) + psi(:,i)**2.0_dp*2.0_dp*energies(i)
+    end do
+     
+    print *, energies(size(energies))
+    dbeta = derivative(bas,beta)
+    f(:) =  - dbeta(:) + energies(size(energies))*dnx(:)
+    !f(:) = !energies(size(energies))*dnx(:)!f(:)
+    do i=1, size(bas%x)
+        do j=1, size(bas%x)-i+1
+            vp(i) = vp(i)- nx(i+j-1)*f(i+j-1)*dx
+        end do
+    end do
+    !vp(:) = 0.0_dp
+    !vp(:) = nx(:)*beta(:)*0.5_dp
+    open(12, file ="vp3.dt")
+    do i=1, size(bas%x)
+        write(12, *) bas%x(i), beta(i)
+    end do
+    
+    close(12)
+    vp(:) = 2.0_dp/(nx(:)+0.001_dp)**2.0_dp * vp(:)
+    do i=1,999
+        vp(i+1001) = vp(1000-i)
+    end do
+        vp(1001) = vp(1000)
+    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 
+    open(12, file ="vp4.dt")
+    do i=1, size(bas%x)
+        write(12, *) bas%x(i), vp(i)
+    end do
+    close(12)
     
+end subroutine
 subroutine get_pauli_beta_exact(H,bas,energies,db,nx)
     real(dp),intent(inout)  :: H(:,:,:)
     real(dp),intent(in)     :: nx(:,:),db,energies(:)

main.f90

@@ -9,7 +9,8 @@ integer :: npeak
 type(eigres_t) :: sol(7), sol2
 type(basis_t) :: testb
 type(dft_t)     :: res, res1(5)
-real(dp):: H(2,100,100),a,b,L,V0,nx(2000),hart(5,2000), pi, psi(2000,2),rhobeta(5,2000),db,dnx(5,2000),oldrbeta(5,2000),cnt,preve
+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):: ergs(5)
 integer :: i,j,k,norbs,nelectrons
 norbs = 5
 nelectrons=10
@@ -18,23 +19,42 @@ pi = 4.0_dp*atan(1.0_dp)
 npeak = 10
 a=1.0_dp
 b=1.0_dp/6.0_dp
-cnt = 0.02_dp
+!cnt = 0.02_dp
 L = 8.0_dp
 V0 = 100.0_dp
 H(:,:,:) = 0.0_dp
-testb = basis_t(nelectrons,10,2000,L)
+testb = basis_t(nelectrons,100,2000,L)
 nx(:) = 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
+    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)
 end do
 
 
-res1(2) = solve_dft(testb,nx,hart(1,:),.false.,norbs,.true.,10)
+!res1(2) = solve_dft(testb,nx,hart(1,:),.false.,norbs,.true.,10)
+
 !print *, sum((res1(1)%nx(:) -res1(2)%nx(:))**2*(testb%x(2)-testb%x(1)) )
-print *, res1(2)%energies
+!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,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)
+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_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)
@@ -44,7 +64,7 @@ call get_pauli_exact(H(1,:,:),testb,res1(2)%psi,res1(2)%energies,res1(2)%nx)
 !print *, res1(2)%energies
 open(12, file="rho.dt")
 do i=1, size(testb%x)
-    write(12, *) testb%x(i),res1(2)%nx(i)
+    write(12, *) testb%x(i),nx(i),nx2(i),nx3(i)
 end do
 close(12)