MODULE calcpotential
  USE GlobalMath
  USE excor
  USE atomdata

  IMPLICIT NONE




CONTAINS
  SUBROUTINE potential(Grid,Pot,ecoul,etxc,eexc)
    !  calculate veff for density functional theory from electron density
    !  nz is nuclear charge
    !  den(n) is electron density * (4*pi*r**2)
    !  rv(n) is returned as veff * r
    !  q is total electron charge
    !  ecoul and etxc are coulomb and exchange-correlation contributions
    !    to the total energy
    !  eexc is the total exchange energy (int(den*exc))
    !  v0=v(0)
    !  v0p=dv/dr(0)
    !  a(n), b(n), c(n), d(n) are work arrays
    TYPE(GridInfo), INTENT(IN) :: Grid
    TYPE(PotentialInfo), INTENT(INOUT) :: Pot
    REAL(8), INTENT(INOUT) :: ecoul,etxc,eexc

    REAL(8), POINTER :: den(:),rv(:)

    INTEGER :: n,nz
    REAL(8) :: h,q,v0,v0p
    REAL(8) :: fpi,r
    REAL(8),  ALLOCATABLE :: tmp(:),tmp1(:)
    INTEGER :: i,j,k

    n=Grid%n
    h=Grid%h
    den=>Pot%den
    rv=>Pot%rv
    q=Pot%q
    nz=Pot%nz

    write(6,*) 'in potential ', nz,q,n,h
    write(6,*) den(1),den(2)
    fpi = 4 * pi      !fpi=16*datan(1.d0)

    ALLOCATE(tmp(n),tmp1(n),stat=i)
    IF (i/=0) THEN
       WRITE(6,*) 'Error in potential allocation ', i,n
       STOP
    ENDIF

    !
    !  Coulomb contribution
    !
    CALL poisson(n,h,q,den,rv,ecoul)
    !
    !  Exchange-correlation contribution  exc-vxc
    !
    CALL exch(n,h,den,tmp,etxc,eexc)
    !
    rv(1:n)=rv(1:n)+tmp(1:n)
    !
    !  add in nuclear contribution and calculate v0 and v0p
    !
    DO i=1,15
       r=i*h
       tmp(i)=rv(i)/r
    ENDDO
    CALL nderiv(h,tmp,tmp1,15,i)
    IF (i/=0) THEN
       WRITE(6,*) 'Error in potential call to nderiv ',i,n
       STOP
    ENDIF
    v0=tmp(3)+3*(tmp(1)-tmp(2))      ;Pot%v0=v0
    v0p=tmp1(3)+3*(tmp1(1)-tmp1(2))  ;Pot%v0p=v0p
    WRITE(6,*) 'v0, v0p = ', v0,v0p
    DO i=1,n
       rv(i)=rv(i)-2*nz
    ENDDO
    DEALLOCATE(tmp,tmp1)
  END SUBROUTINE potential


  SUBROUTINE poisson(n,h,q,den,rv,ecoul)
    !  use Numerov algorithm to solve poisson equation
    !  den(n) is electron density * (4*pi*r**2)
    !  rv(n) is returned as electrostatic potential * r
    !  q is total electron charge
    !  a(n) and b(n) are work arrays
    !  ecoul is the coulomb interation energy

    INTEGER, INTENT(IN)::n
    REAL(8), INTENT(IN):: h,q,den(:)
    REAL(8), INTENT(INOUT) :: rv(:),ecoul

    REAL(8), ALLOCATABLE :: a(:),b(:)
    REAL(8) :: sd,sl,ss2,th
    INTEGER :: i

    ALLOCATE(a(n),b(n),stat=i)
    IF (i/=0) THEN
       WRITE(6,*) 'Error in allocating arrays in poisson',i,n
       STOP
    ENDIF

    rv=0
    sd=2
    sl=-1
    DO i=1,n
       a(i)=h*den(i)/(6*i)
    ENDDO
    rv(1)=10*a(1)+a(2)
    DO i=2,n-1
       rv(i)=10*a(i)+a(i+1)+a(i-1)
    ENDDO
    rv(n)=10*a(n)+a(n-1)+2*q
    ! use Thomas's algorithm for inverting matrix
    !   Dale U. von Rosenberg, "Methods for the Numerical Solution of
    !     Partial Differential Equations, Am. Elsevier Pub., 1969, pg. 113
    a(1)=sd
    ss2=sl*sl
    DO i=2,n
       a(i)=sd-ss2/a(i-1)
    ENDDO
    b(1)=rv(1)/sd
    DO i=2,n
       b(i)=(rv(i)-sl*b(i-1))/a(i)
    ENDDO
    rv(n)=b(n)
    DO i=n-1,1,-1
       rv(i)=b(i)-sl*rv(i+1)/a(i)
    ENDDO
    !
    !
    !  calculate ecoul
    !
    DO i=1,n
       a(i)=den(i)*rv(i)/i
    ENDDO
    th=1
    ecoul=0.5d0*overint(n,th,a)
    WRITE(6,*) 'ecoul = ',ecoul
    DEALLOCATE(a,b)
    !
    RETURN
  END SUBROUTINE poisson

END MODULE calcpotential