! Simple nuclear Hartree-Fock program. 
      program simplehf
      implicit none
c
      integer npoints
      real t0,t3
      parameter (npoints = 40)
      parameter (t0 = -1090.0, t3 = 17488.0)
      integer  i, j, maxiter
      real dr, fourpi, hbrm, hob, xdamp
      real  spe, ke, speold, sum
      real wavefunction(npoints),density(npoints),potential(npoints)
      real r_grid(0:npoints) , scr1(npoints)
       maxiter = 20
       sum=0.0
       dr = 0.2
       fourpi = 12.56637061432
       hbrm = 41.437
       hob = 1.40
       xdamp = 0.0001 ! damping parameter
c
c
c     Initialize r-grid
      do i=0,npoints-1
         r_grid(i)=real(i)*dr
      end do
      do i=1,npoints
           scr1(i)=0.0
      end do

c     Initialize the wave function with a gaussian
      do i=1,npoints
         wavefunction(i)=(r_grid(i-1)/hob)
     &  *exp(-(r_grid(i-1)/hob)**2/2.0)
         sum=sum+wavefunction(i)*wavefunction(i) !sum of square
      end do

c     normalize wavefunction
      do i=1,npoints
         wavefunction(i) = wavefunction(i)/sqrt(dr*sum)
      end do

c     Start HF iteration here
      speold = 1.0e10 ! large positive number
      do i=1,maxiter
calculate density
         do j=2,npoints
            density(j)=4.0*(wavefunction(j)/r_grid(j-1))**2/fourpi
         end do
c special at the origin based on L'Hopital's rule
       density(1) =4.0*
     & ((4.0*wavefunction(2)/3.0-wavefunction(3)/6.0)/dr)**2/fourpi
        do j=1,npoints
           potential(j)=3.0*t0/4.0*density(j)+3.0*t3/16.0*density(j)**2
        end do
c act with kinetic energy operator on wavefunction using 3-pt. formula
        scr1(1) = hbrm*wavefunction(1)/dr**2
        do j=2,npoints-1
           scr1(j) = -hbrm/2.*(wavefunction(j-1)-2.*wavefunction(j)+
     &     wavefunction(j+1))/dr**2
        end do
        scr1(npoints) = -hbrm/2*(wavefunction(npoints-1)-
     &  wavefunction(npoints))/dr**2
        sum=0.0
        do j=1,npoints
           sum=sum+scr1(j)*wavefunction(j)
        end do
        ke=dr*sum
        do j=1,npoints
           scr1(j)=scr1(j)+potential(j)*wavefunction(j)
        end do
        sum=0.0
        do j=1,npoints
           sum=sum+scr1(j)*wavefunction(j)
        end do
        spe=dr*sum   
c        ke=dr*sum 
c  print out some information every 20 iterations
        if (mod(i,20).EQ.0) then
           write(*,*) "It,ke,spe,conv=",i,ke,spe,abs(spe-speold)
        end if     
        speold=spe
c  update wavefunction
        do j=1,npoints
           wavefunction(j)=wavefunction(j)
     &     -xdamp*(scr1(j)-spe*wavefunction(j))
        end do
        sum=0.0
        do j=1,npoints
           sum=sum+wavefunction(j)*wavefunction(j)
        end do
        do j=1,npoints
           wavefunction(j)=wavefunction(j)/sqrt(dr*sum)
        end do
      end do
  ! finally, print the density and potential
      do i=1,npoints
         write(15,*) r_grid(i-1),wavefunction(i),density(i),
     &   fourpi*r_grid(i-1)**2*density(i)
         write(16,*) r_grid(i-1),potential(i)
      end do
      end ! program simplehf

