PROGRAM Hartree
IMPLICIT NONE
! Calculates the ground state energy of two electron atoms
! See "Physics of Atoms and Molecules" by B.H. Brandsen and C.J. Joachain
INTEGER N, NMAX, ITER, KITER
! N is the radial distance counter
! NMAX+1 is the number of radial points starting a r=0 and ending at r=RMAX
! ITER is the maximum number of interations allowed
! KITER is the iteration counter
REAL(KIND=8) EPS, ZETA, W, DELTA, EAV, EK, EN, Z, RMAX, TOL, DEPS, ESS, H, PNMAX, SUM
! EPS is the orbital energy
! ZETA is the number of protrons in the nucleus
! W is Numer of electrons
! DELTA is the error estimate
! EAV is the HF energy
! EK is the kinetic energy
! EN is the nuclear energy term
! Z is the number of protrans in the nucleus
! RMAX is the maximum radial distance in units of Bohr
! TOL is the desired error tolerance for the energy (Hartree)
! DEPS is the current error estimate
REAL(KIND=8), ALLOCATABLE:: P(:),V(:)
! P(0:NMAX) is an array containing the radial wavefunction
! V(0:NMAX) is an array containing the radial potential
WRITE (6,*) ' INPUT: NMAX,EPS,ZETA,W,Z,RMAX,TOL,ITER,DEPS'
READ (5,*) NMAX, EPS, ZETA , W, Z, RMAX, TOL, ITER, DEPS

ALLOCATE(P(0:NMAX), V(0:NMAX))
H=RMAX/NMAX
DO N=0,NMAX
    P(N)=2.0D0*ZETA**1.5D0*DEXP(-ZETA*N*H)*N*H
ENDDO

WRITE (6,'(5A22)') "EAV","EPS","EK","EN","ESS"
DO KITER=1,ITER
    DELTA=DEPS
    SUM=1.D0
    V(NMAX)=2.D0*(W-1.D0-Z)/RMAX
    DO N=NMAX-1,1,-1
        SUM=SUM-P(N+1)**2*H
        V(N)=V(N+1)+2.D0/(N*(N+1))*(SUM*(W-1.D0)-Z)/H
    ENDDO
300 SUM=P(1)**2*H/2.0D+00
    DO N=1,NMAX-1
        P(N+1)=(2.0D0*P(N)-P(N-1)+H*H/12.0D0*(10.0D0*(V(N)-EPS)*P(N)+ (V(N-1)-EPS)*P(N-1)))/(1.0D0-(V(N+1)-EPS)*H*H/12.0D0)
        SUM=SUM+P(N+1)**2*H
    ENDDO
    IF (PNMAX*P(NMAX).LT.0.D0) DELTA=DELTA/2.0D+00
    PNMAX=P(NMAX)
    EPS=EPS+DSIGN(1.D0,PNMAX)*DELTA
    IF (DELTA.GT.TOL) GO TO 300
    EN=0.D0
    ESS=0.D0
    DO N=1,NMAX
        P(N)=P(N)/DSQRT(SUM)
        EN=EN-2.D0*Z/(N*H)*P(N)**2*H
        ESS=ESS+(2.D0*Z/(N*H)+V(N))*P(N)**2*H
    ENDDO
    EAV=W*EPS-ESS*(W-1.D0)
    EK=EPS-EN-ESS*(W-1.D0)
    WRITE (6,'(5ES22.14)') EAV,EPS,EK,EN,ESS
ENDDO

DEALLOCATE(P,V)
STOP 'Execution of Hartree is complete'
ENDPROGRAM Hartree