joitain korjauksia

energyminimizer.f90

@@ -7,7 +7,6 @@ module energyminimizer
 	use derivatives, only: rDiff3, rDiff1, rDiff2, Weizsacker
 	use hartreepot, only: hartree
 	use xc, only: LDAvX,LDAvC
-	use kedf_nagy, only: CalculateNagyPotential
 	implicit none
 	public tsolver, ConjugateGradientMethod
 

kedf_nagy.f90

@@ -246,8 +246,8 @@ contains
 		!vp(:) = 1.0_dp/6.2_dp * vp(:)
 
 		!vp = 4*3.141592653 * r(:)**2 * vp(:)
-		!print *, trapz7(rho(:) * (vp(:) + Weizsacker(rho)) - orbenergies(NumberOfOrbs)*rho(:)&
-		!+(orbenergies(1)*2*orbdens(:,1)+orbenergies(2)*2*orbdens(:,2)),1,Nmesh)
+		print *, trapz7(rho(:) * (vp(:) + Weizsacker(rho)) - orbenergies(NumberOfOrbs)*rho(:)&
+		+(orbenergies(1)*2*orbdens(:,1)+orbenergies(2)*2*orbdens(:,2)),1,Nmesh)
 		deallocate(qinputs)
 		deallocate(qdotinputs)
 		deallocate(qinputsPrev)

kineticenergies.f90

@@ -23,8 +23,8 @@ contains
 		allocate(occupations(NumberOfOrbs))
 		allocate(l(NumberOfOrbs))
 		allocate(orbenergies(NumberOfOrbs))
-		Z=9
-		open(unit=11,file="results_Fe_5001_fullorbs")
+		Z=26
+		open(unit=11,file="results_F")
 		orbdens = 0.0_dp
 		do i=1, Nmesh
 			read(11,*) grid(i), rho(i), orbdens(i,1), orbdens(i,2), orbdens(i,3),orbdens(i,4),orbdens(i,5),orbdens(i,6), &
@@ -48,9 +48,9 @@ contains
 		!orbdens(:,8) =  (1.0_dp/trapz7(orbdens(:,8),1,Nmesh)) * orbdens(:,8)
 		!orbdens(:,9) =  (1.0_dp/trapz7(orbdens(:,9),1,Nmesh)) * orbdens(:,9)
 		vks(:) = hartree(rho) + LDAvC(rho) + LDAvX(rho) - 26.0_dp/grid(:)
-		print *, NagyKineticEnergy(rho(:), orbdens(:,:), vks(:), occupations, l)
+		!print *, NagyKineticEnergy(rho(:), orbdens(:,:), vks(:), occupations, l)
 		!print *, trapz7(Weizsacker(rho(:)) *rho, 1, Nmesh)
-		!call PauliPotential(initrho, rho, orbdens, orbenergies, occupations, l )
+		call PauliPotential(initrho, rho, orbdens, orbenergies, occupations, l )
 	end subroutine TestDerivatives
 
 

main_Be.f90

@@ -49,7 +49,7 @@ program main_Be
   aMesh=2.7e6_dp     ! a = r(last) / r(first): controls the change in point density
   rmin=1.0e-7_dp ! Location of the first point -7 stand
   rmax=50.0_dp   ! Location of the last point
-  Nmesh=5001     ! Number of points
+  Nmesh=2001    ! Number of points
   !
   ! Accuracy of the finite difference approximations:
   ! stencil = number of points in the finite differences
@@ -80,10 +80,10 @@ program main_Be
   vXC=LDAvX(RhoOld)+LDAvC(RhoOld)
   vKS=vExt+vH+vXC
   !
-  ConvTarget=1.0e-5
+  ConvTarget=1.0e-8
   allocate(initrho(Nmesh))
   ! T/he main iteration loop can now begin
-  goto 10
+ ! goto 10
   do iterMain=1,itersMainMax
      print *,''
      print *,'Iteration number ',iterMain
@@ -174,7 +174,7 @@ program main_Be
      write(*,*) ''
   end do
   !     xclb      = 2.871234 * trapz7(RhoNew(:,1) **(5/3),1,Nmesh)
-  10 CONTINUE
+  !!10 CONTINUE
   allocate(testvp(Nmesh))
 
   !testvp = vpsolver(rhoNew(:,1), orbs(1,:), orbs(2,:), e1s, muNew)
@@ -183,27 +183,27 @@ program main_Be
   ! Save calculated quantities to disk.
   filename = "./results"
   allocate(outArray(Nmesh,4))
-  allocate(vwp(Nmesh))
-  allocate(muar(Nmesh))
-  allocate(betaterm2(Nmesh))
-  allocate(dTdn(Nmesh))
-  betaterm2 = 1.0_dp
-  initrho = RhoNew(:,1)!1.0_dp/50.0_dp!4*3.141592*r(:)**3*exp(-2*(r(:) + 5.0_dp))!(8.0_dp*exp(-2.0_dp*R*Z)*R**2.0_dp*Z**3.0_dp+&
+  !allocate(vwp(Nmesh))
+  !allocate(muar(Nmesh))
+  !allocate(betaterm2(Nmesh))
+  !allocate(dTdn(Nmesh))
+  !betaterm2 = 1.0_dp
+  !initrho = RhoNew(:,1)!1.0_dp/50.0_dp!4*3.141592*r(:)**3*exp(-2*(r(:) + 5.0_dp))!(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)!4.0_dp/trapz7(betaterm2,1,Nmesh)4
 
 !  initrho = initrho*1.0_dp/trapz7(initrho,1,Nmesh)
-  dTdN = Weizsacker(RhoNew(:,1)) + vPNew(:,1)
-  betaterm2 = initrho
+  !dTdN = Weizsacker(RhoNew(:,1)) + vPNew(:,1)
+  !betaterm2 = initrho
   !betaterm2 = SlaterDensity(8.0_dp, 5.3_dp)
   !call CalculateNagyPotential(RhoNew(:,1), vks(:), betaterm2(:))
   !call ChemicalPotential(vks(:) + betaterm2(:),sqrt(RhoNew(:,1)), 4, muNew)
   !print *, muNew
   !print *, NagyKineticEnergy(RhoNew(:,1), vks(:))
-  vwp = Weizsacker(RhoNew(:,1))
-  muar = muNew
+  !vwp = Weizsacker(RhoNew(:,1))
+  !muar = muNew
   !print *, trapz7(RhoNew(:,1) * dTdn, 1, Nmesh), trapz7(betaterm2(:) * RhoNew(:,1), 1, Nmesh)
   outArray(:,1)=r
-  call TestDerivatives(initrho,betaterm2)
+  !call TestDerivatives(initrho,betaterm2)
  ! call DirectGradientMethod(initrho(:))
   !call MinimizeSlaterGrid(RhoNew(:,1),vks(:))
   !call CalculateTotalEnergy(energy, betaterm2)
@@ -216,8 +216,8 @@ program main_Be
 !print *, muNew
   !outArray(:,3)  Psol(:,1)
   !outArray(:,2)=vPNew(:,2)
-  outArray(:,2) = initrho(:)
-  outArray(:,3)= betaterm2(:)
+  !outArray(:,2) = Qdot(:,1)
+  !outArray(:,3)= betaterm2(:)
   !outArray(:,4)= orbs(:,2)
   !outArray(:,2) = vks(:)
     !outArray(:,5)=

minimize_directgradient.f90

@@ -33,15 +33,15 @@ contains
 		l(:) = 0
 		rhos(:) = sqrt(rho(:))
 		rhoOld(:) = rhos(:)
-		open(unit=11,file="results_Be_2001_fullorb")
-		do i=1, Nmesh
-			read(11,*) dummy(i),realrho(i),orb1(i),orb2(i)
-		end do
-		close(11)
-		orbs(:,1) = (1.0_dp/trapz7(orb1(:),1,Nmesh)) * orb1(:)
-		orbs(:,2) = (1.0_dp/trapz7(orb2(:),1,Nmesh)) * orb2(:)
-
-		print *, trapz7(rhoOld, 1, Nmesh)
+		!open(unit=11,file="results_Be_2001_fullorb")
+		!do i=1, Nmesh
+		!	read(11,*) dummy(i),realrho(i),orb1(i),orb2(i)
+		!end do
+		!close(11)
+		!orbs(:,1) = (1.0_dp/trapz7(orb1(:),1,Nmesh)) * orb1(:)
+		!orbs(:,2) = (1.0_dp/trapz7(orb2(:),1,Nmesh)) * orb2(:)
+
+		print *, trapz7(rhoOld, 1, Nmesh), "hey"
 		do i=1, 10000
 		
 			call CalculatePotential(rhos**2, pot, orbs, occupations, l, realrho)
@@ -112,7 +112,7 @@ contains
 		!call CalculateNagyPotential(rho,vks,T)
 		!call PauliPotential(vp, realrho, orbs,orbenergies, occupations, l )
 
-		T(:) =  (1/r(:)**2)* rDiff1(r(:)**2*rDiff1(sqrt(realrho)))+ 2*sqrt(rho)*vp(:)
+		T(:) =  (1/r(:)**2)* rDiff1(r(:)**2*rDiff1(sqrt(realrho)))!+ 2*sqrt(rho)*vp(:)
 		if(any(isnan(T))) then
 			print *, "weizsäcker nan"
 		end if

solverhovp.f90

@@ -95,9 +95,10 @@ contains
     e1sTarget=1.0e-6_dp
     e1sShift=1.0e-5_dp
     e1sMultiplier=1.0_dp
-
+	e1s = -3.856411_dp
     call AdamsOutwardCTP(Pout,Pdotout,Pgdotout,Qout,Qdotout,Qgdotout,orbs,vPout,mu,vKS,ctp,e1s,gamma,start)
         !call AdamsInward(Pin,Pdotin,Qin,Qdotin,vPin,mu,vKS,ctp)
+    goto 10
     vPmin=vPout(ctp,1)
     if (prnt==1) then
        print *,'e1s,vPmin = ',e1s,vPmin
@@ -170,6 +171,7 @@ contains
     ! 
     ! Now we are able to find the correct mu
     ! 
+    10 CONTINUE
     do iterAdams=1,iterAdamsMax
        !print *,'Adams iteration number ',iterAdams
        call AdamsOutward(Pout,Pdotout,Pgdotout,Qout,Qdotout,Qgdotout,orbs,vPout,mu,vKS,ctp,e1s,gamma,start)