## A Simple Hartree-Fock program

Send questions to Paul Stevenson
```
! 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

```