A Simple Hartree-Fock program

! Simple nuclear Hartree-Fock program.  Paul Stevenson 2003.
! http://www.ph.surrey.ac.uk/~phs3ps/simple-hf.html

program simplehf

  implicit none

  integer :: i, j

  integer, parameter :: maxiter = 1000
  integer, parameter :: npoints = 40
  real, parameter :: dr = 0.2
  real, parameter :: fourpi = 12.56637061432
  real, parameter :: hbrm = 41.437
  real, parameter :: hob = 1.4

  ! damping parameter
  real, parameter :: xdamp = 0.0001

  ! storage of single particle state
  real, dimension(npoints) :: wavefunction

  ! density
  real, dimension(npoints) :: density

  ! HF potential 
  real, parameter :: t0 = -1090.0, t3 = 17488.0
  real, dimension(npoints) :: potential

  ! a vector with the r-grid points
  real, dimension(npoints) :: r_grid 

  ! some workspace
  real, dimension(npoints) :: scr1

  ! single particle potential and kinetic energies
  real :: spe, ke, speold

  ! Initialize r-grid
  r_grid = (/ (i*dr,i=0,npoints-1) /)

  ! Initialize the wave function with a gaussian
  wavefunction = (r_grid/hob) * exp(-(r_grid/hob)**2/2)

  ! normalize wavefunction
  wavefunction = wavefunction / sqrt(dr*sum(wavefunction**2))

  ! Start HF iteration here
  speold=HUGE(1.0)
  do i=1,maxiter
     ! calculate density
     density(2:npoints) = 4*( wavefunction(2:) / r_grid(2:) )**2 / fourpi
     ! special at the origin based on L'Hôpital's rule
     density(1) = 4*( (4*wavefunction(2)/3-wavefunction(3)/6) /dr)**2/fourpi

     potential = 3*t0/4 * density + 3*t3/16 * density**2

     ! 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

     ke = dr * sum(scr1*wavefunction)

     scr1 = scr1 + potential*wavefunction
     spe = dr * sum(scr1*wavefunction)
     
     ! print out some information every 20 iterations
     if(mod(i,20)==0) write(6,*) "It,ke,spe,conv=",i,ke,spe,abs(spe-speold)
     
     speold=spe
     ! update wavefunction
     wavefunction = wavefunction - xdamp * (scr1-spe*wavefunction)
     wavefunction = wavefunction / sqrt(dr*sum(wavefunction**2))

  end do
  
  ! finally, print the density and potential
  do i=1,npoints
     write(15,*) r_grid(i),density(i),fourpi*r_grid(i)**2*density(i)
     write(16,*) r_grid(i),potential(i)
  end do

end program simplehf