diff --git a/analyze_vp_error.py b/analyze_vp_error.py new file mode 100644 index 0000000000000000000000000000000000000000..4df0dfbb416dad428600687821cb9a3abc3f8161 --- /dev/null +++ b/analyze_vp_error.py @@ -0,0 +1,15 @@ +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 diff --git a/hamilton.f90 b/hamilton.f90 index dc46a55a4239c87ad0af61676cd24943e148c977..953b5c46a176d71c8355318a8201f049c9b5eedc 100644 --- a/hamilton.f90 +++ b/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(:) diff --git a/main.f90 b/main.f90 index 9384f5b70d359afaee0b9411af9849f6d4c56f33..bf70ef427b823f2a69d599baa0f348080e34045a 100644 --- a/main.f90 +++ b/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)