c---------------------------------------------------------------
      program atom_slaterx
c---------------------------------------------------------------
      implicit none
      integer mshx,nwfx
      parameter (mshx = 3001, nwfx=10)
      real*8 pi, fpi
      parameter (pi=3.14159265358979d0, fpi=4.d0*pi)

      integer mesh
      integer n, l, i, iter
      real*8 zmesh, rmax, xmin, zeta, dx, de, tr2, beta, tol
      real*8 ehar, extot, vxtot, evion, ekin, etot, etot1, sum_eigenv
      real*8 r(0:mshx), sqr(0:mshx), r2(0:mshx)
      real*8 yy(0:mshx), vpot(0:mshx), work(0:mshx), ee
      real*8 rho(0:mshx), vscf(0:mshx), vhx(0:mshx), deltav(0:mshx)
      integer nwfc, nn(nwfx), ll(nwfx)
      real*8 e(nwfx), occ(nwfx), y(0:mshx,nwfx)
c
c read atomic charge, and other parameters
c
      write(*,*) 
      write(*,*) 'Hartree+SlaterX calculation for generic atom'
      write(*,*) 
      write (*,'(a,$)') ' ATOMIC CHARGE zeta [ >=2 ] '
      read(*,*) zeta
      write(*,*) zeta
      if (zeta.lt.2.0) stop 'zeta should be .gt. 2'
      write (*,'(a,$)') ' MIXING PARAMETER beta [0.0-1.0] '
      read(*,*) beta
      write(*,*) beta
      if (beta.lt.0.0 .or. beta.gt.1.0) stop 'beta not in [0.0-1.0]'
      write (*,'(a,$)') ' SCF ACCURACY     = ', tol
      read(*,*) tol
      write(*,*) tol
      if (tol.lt.0.0) stop 'tol shold be positive'
c
c initialize logarithmic mesh
c
      zmesh=  zeta 
      rmax = 100.d0
      xmin = -6.d0
      dx   =  0.01d0
      call do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
c
c initialize the potential
c
      call init_pot(zeta,r,mesh,vpot)
c
c read the configuration for the atom
c
      write(*,*) ' nwfc >> '
      read(*,*) nwfc
      if (nwfc.gt.nwfx) stop 'nwfc.gt.nwfx'

      do i=1,nwfc
         write(*,*) 'nn,ll,occ >>'
         read(*,*) nn(i),ll(i),occ(i)
         if (nn(i).lt.ll(i)+1) then
            write(*,*) 'error in main: n.lt.l -> wrong number of nodes '
            stop
         end if
      end do

      do i=0,mesh
         vscf(i) = vpot(i)
      end do

      iter = 0
    1 continue ! this is the entry point of the SCF cycle
      iter = iter + 1
      write(*,*) '####################################################'
      write(*,'(a,i3)') ' SCF iteration # ',iter
c
c solve the schroedinger equation in radial coordinates by Numerov method
c
      do i=1,nwfc
         call solve_sheq(nn(i),ll(i),e(i),mesh,dx,r,sqr,r2,
     +                   vscf,zeta,y(0,i),work)
      end do
c
c calculate the charge density from the wfc
c
      call rho_of_r(mesh,mshx,r,r2,y,occ,nwfc,rho)
c
c calculate the Hartree + Exchange potential and energy
c
      call v_of_rho(mesh,mshx,dx,r,r2,rho,vhx,ehar,extot,vxtot)
c
c calculate the kinetic energy and the energy in the external potential
c 
      evion = 0.d0
      sum_eigenv = 0.d0
      do i=1,nwfc
         sum_eigenv = sum_eigenv + occ(i) * e(i)
      end do
      ekin = sum_eigenv
      do i=0,mesh
         evion = evion + vpot(i) * rho(i) * fpi * r2(i) * r(i) * dx
         ekin  = ekin - vscf(i) * rho(i) * fpi * r2(i) * r(i) * dx
      end do

      de = 0.d0
      tr2 = 0.d0
      do i=0,mesh
         deltav(i) = vpot(i) + vhx(i) - vscf(i) 
         vscf(i) = vscf(i) + beta * deltav(i)
         de = de + deltav(i) * rho(i) * fpi * r2(i) * r(i) * dx
         tr2 = tr2 + deltav(i)*deltav(i)*r(i)*dx
      end do
c
c write out the eigenvalue energy to be compared with the external potential
c
      etot = sum_eigenv - ehar + extot - vxtot + de
      etot1= ekin + evion + ehar + extot

      do i=1,nwfc
         write (*,'(a,2i3,a,f15.6)') 'eigenvalue ', nn(i),ll(i),'  = ',
     +                                              e(i)
      end do
      write (*,'(a,f15.6)') 'Eigenvalue energy       ', sum_eigenv
      write (*,'(a,f15.6)') 'Kinetic energy          ', ekin
      write (*,'(a,f15.6)') 'External pot. energy    ', evion
      write (*,'(a,f15.6)') 'Hartree energy          ', ehar
      write (*,'(a,f15.6)') 'Exchange energy         ', extot - vxtot
      write (*,'(a,f15.6)') 'Variational corr. energy', de
      write (*,'(a,2f15.6)')'Total energy            ', etot, etot1
      write (*,'(a,f15.6)') 'Virial check            ',
     +                                        -(evion+ehar+extot)/ekin
      write (*,'(a,1pe15.1)') 'tr2                    ', tr2

      if(tr2.gt.tol) go to 1

      write (*,*) 
      write (*,*) ' SCF Convergence has been achieved'
      write (*,*) 
      write (*,*) ' compute additional single-particle states'
      write (*,*) 
c
c write on file pot.out
c
      open (7,file='pot.out',status='unknown',form='formatted')
      do i =0,mesh
         write (7,*) r(i),vpot(i),vhx(i),vscf(i)
      end do
      close(7)
c
c open wfc file
c
      open (7,file='wfc.out',status='unknown',form='formatted')

    2 continue     
c
c read principal and angular quantum numbers
c
      write (*,'(a,$)') ' n, l >'    
      read(*,*) n,l
      if (n.lt.l+1) then
         write(*,*) 'error in main: n.lt.l -> wrong number of nodes '
         go to 2
      end if
      call solve_sheq(n,l,ee,mesh,dx,r,sqr,r2,vscf,zeta,yy,work)
      write (*,'(a,2i3,a,f15.6)') 'eigenvalue ',n,l,'  = ',ee
      do i=0,mesh
         write (7,*) r(i),yy(i)/sqr(i), yy(i)*sqr(i), ee
      end do
      write (7,*) 0.d0,0.d0, 0.d0, ee

      go to 2

      end 
c---------------------------------------------------------------------
      subroutine solve_sheq(n,l,e,mesh,dx,r,sqr,r2,vpot,zeta,y,f)
c---------------------------------------------------------------------
c
c solve the schroedinger equation in radial coordinates by Numerov method
c
      implicit none
      integer mesh, n,l,j
      integer i, icl, nodes, ncross, kkk
      real*8 dx, ddx12, sqlhf, x2l2, zeta
      real*8 r(0:mesh), sqr(0:mesh), r2(0:mesh), 
     +       vpot(0:mesh), y(0:mesh), f(0:mesh)
      real*8 ycusp, dfcusp, fac, norm, eup, elw, e, de

      ddx12=dx*dx/12.d0
      sqlhf = (l+0.5d0)**2
      x2l2 = 2*l+2
c
c set initial lower and upper bounds to the eigenvalue
c
      eup=vpot(mesh)
      elw=eup
      do  i=0,mesh
          elw=min(elw,sqlhf/r2(i)+vpot(i))
      end do
      if (eup-elw.lt.1.d-10) then
         write (*,*) elw, eup
         write (*,*) 'solve_sheq: lower and upper bounds are equal'
         stop
      end if
      e = 0.5d0 * (elw + eup)

      kkk = 0
  1   continue ! this is the entry point for the solution at fixed energy
      kkk = kkk + 1
c
c set up the f-function and determine the position of its last
c change of sign
c f < 0 (approximatively) means classically allowed   region
c f > 0         "           "        "      forbidden   "
c
      icl = 2
      f(0) = ddx12 *( sqlhf + r2(0) * (vpot(0)-e) )
      do i=1,mesh
         f(i) = ddx12 * ( sqlhf + r2(i) *(vpot(i)-e) )
         if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i
      end do

      if(icl.ge.mesh-2) then
c         write (*,*) icl, mesh
c         write (*,*) 'error in solve_sheq: last change of sign too far'
         eup = e
         e = 0.5d0 * (eup+elw)
         go to 1
      end if

      do i=0,mesh
         f(i)=1.0d0-f(i)
         y(i)=0.d0
      end do 
c
c determination of the wave-function in the first two points 
c
      nodes = n - l - 1
      y(0) = r(0)**(l+1) *(1.d0 - 2.d0*zeta*r(0)/x2l2) / sqr(0)
      y(1) = r(1)**(l+1) *(1.d0 - 2.d0*zeta*r(1)/x2l2) / sqr(1)
c
c outward integration and count number of crossings
c
      ncross=0
      do i =1, icl-1
         y(i+1)=((12.d0-10.d0*f(i))*y(i)-f(i-1)*y(i-1))/f(i+1)
         if ( y(i) .ne. sign(y(i),y(i+1)) ) ncross=ncross+1
      end do
      fac = y(icl) 
c
c check number of crossings
c
      if (ncross.ne.nodes) then
         if (kkk.gt.100) then
            write (*,*) kkk, e, elw, eup, ncross, nodes, icl
            write (*,*) ' error in solve_sheq: too many iterations'
            stop
         end if
         if (ncross.gt.nodes) then
            eup = e
         else 
            elw = e
         end if
         e = 0.5d0 * (eup+elw)
         go to 1
      end if
c
c determination of the wave-function in the last two points 
c assuming y(mesh+1) = 0 and y(mesh) = dx
c
      y(mesh) = dx
      y(mesh-1) = (12.d0-10.d0*f(mesh))*y(mesh)/f(mesh-1) 
c
c inward integration 
c
      do i = mesh-1,icl+1,-1
         y(i-1)=((12.d0-10.d0*f(i))*y(i)-f(i+1)*y(i+1))/f(i-1)
         if (y(i-1).gt.1.d10) then
            do j=mesh,i-1,-1
               y(j) = y(j)/y(i-1)
            end do
         end if
      end do
c
c rescale function to match at the classical turning point (icl)
c
      fac = fac/y(icl)
      do i=icl,mesh
         y(i) = y(i)*fac
      end do
c
c normalize on the segment
c
      norm = 0.d0
      do i=1,mesh
         norm = norm + y(i)*y(i) * r2(i) * dx
      end do
      norm = sqrt(norm)
      do i=0,mesh
         y(i) = y(i) / norm
      end do
c
c find the value of the cusp at the matching point (icl)
c
      i = icl
      ycusp = (y(i-1)*f(i-1)+f(i+1)*y(i+1)+10.d0*f(i)*y(i)) / 12.d0
      dfcusp = f(i)*(y(i)/ycusp - 1.d0)
c
c eigenvalue update using perturbation theory
c
      de = dfcusp/ddx12 * ycusp*ycusp * dx 
      if (de.gt.0.d0) elw = e
      if (de.lt.0.d0) eup = e
      e = max( min (e+de,eup),elw)
c
c the following update is possibly safer
c
c      if (abs(de).gt.eup-elw) then
c         e = 0.5d0*(eup+elw) 
c      else 
c         e = e + de
c      end if
c


      if (kkk.gt.100) then
         write (*,*) kkk, e, de
         write (*,*) ' error in solve_sheq: too many iterations'
         stop
      end if

      if (abs(de).gt.1.d-10) go to 1

c ---- convergenge has been acheived -----
c      write (*,*) ' convergence acheived at iter #',kkk,
c     ,            ' de = ', de

      return
      end
c--------------------------------------------------------------------
      subroutine do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
c--------------------------------------------------------------------
c
c initialize mesh
c
      implicit none
      integer mshx, mesh
      real*8 zmesh,xmin,dx,rmax,r(0:mshx), sqr(0:mshx), r2(0:mshx)

      integer i
      real*8 x

      mesh = (log(zmesh*rmax) - xmin)/dx
      if (mesh.gt.mshx) then
         write (*,*) mesh, mshx
         write (*,*) 'error in routine do_mesh: mesh.gt.mshx '
         stop
      end if

      do i=0,mesh
         x = xmin + dx * i
         r(i)  = exp(x)/zmesh
         sqr(i)= sqrt(r(i))
         r2(i) = r(i) * r(i)
      end do
      write(*,*) ' radial mesh information:'
      write(*,'(a,f10.6,a,f10.6,a,f10.6)') 
     +          '  dx =',dx,', xmin =', xmin,', zmesh =',zmesh
      write(*,'(a,i6,a,f10.6,a,f10.6)')
     +          '  mesh =',mesh,', r(0) =',r(0),', r(mesh) =', r(mesh)
      write(*,*) 
      
      return
      end
c--------------------------------------------------------------------
      subroutine init_pot(zeta,r,mesh,vpot)
c--------------------------------------------------------------------
c
c initialize potential
c
      implicit none

      integer mesh
      real*8 zeta, r(0:mesh), vpot(0:mesh)
      integer i

      open (7,file='pot.out',status='unknown',form='formatted')
      do i =0,mesh
         vpot (i) = - 2.d0*zeta/r(i)
         write (7,*) r(i),vpot(i)
      end do
      close(7)

      return
      end
c--------------------------------------------------------------------
      subroutine rho_of_r(mesh,mshx,r,r2,y,occ,nwfc,rho)
c--------------------------------------------------------------------
c
c compute the charge density of an He-like atom
c
      implicit none
      real*8 pi, fpi
      parameter (pi=3.14159265358979d0, fpi=4.d0*pi)

      integer mesh, mshx, nwfc
      real*8 r(0:mshx), r2(0:mshx), y(0:mshx,nwfc), occ(nwfc)
      real*8 rho(0:mshx)

      integer i,iwf

      do i=0,mesh
         rho(i) = 0.d0
         do iwf=1,nwfc
            rho(i) = rho(i) + 
     +               occ(iwf) * y(i,iwf)**2 * r(i) / (fpi*r2(i))
         end do
      end do

      return
      end 
c--------------------------------------------------------------------
      subroutine v_of_rho(mesh,mshx,dx,r,r2,rho,vhx,ehar,extot,vxtot)
c--------------------------------------------------------------------
c
c compute the Hartree + Exchange potential for an He-like atom
c Exchange effects are computed within the Slater approximation
c
      implicit none
      real*8 e2, pi, fpi
      parameter (e2=2.d0, pi=3.14159265358979d0, fpi=4.d0*pi)

      integer mesh,mshx
      real*8 dx, r(0:mshx), r2(0:mshx), rho(0:mshx)
      real*8 vhx(0:mshx), ehar, extot, vxtot, ex, vx

      integer i
      real*8 charge

c
c calculate the Hartree potential and energy by integrating the 
c electric field generated by the electronic charge. 
c This is done in 2 steps
c
c 1) calculate the charge inside a sphere and fill vhx with the 
c    electric field generated by this charge
c
      charge = 0.d0
      do i=0,mesh
         charge = charge + rho(i) * fpi * r2(i) * r(i) * dx
         vhx(i) = e2*charge/r2(i)
      end do
c
c   (the total cherge is writen in output as a check)
c 
      write (*,'(a,f15.6)') 'Total charge = ', charge
c
c 2) integrate the electric field from +\infty to r to get the potential
c    and integrate V_{Hartree}*rho to get the energy
c
      ehar = 0.d0
      vhx(mesh) = e2 * charge / r(mesh)
      do i = mesh-1, 0, -1
         vhx(i) = vhx(i+1) + vhx(i) * r(i) * dx
         ehar = ehar + vhx(i) * rho(i) * fpi * r2(i) * r(i) * dx
      end do
      ehar = ehar/2.d0
c
c Exchange effects are computed within the Slater approximation
c
      extot = 0.d0
      vxtot = 0.d0
      do i=0,mesh
         call slaterx(rho(i),ex,vx)
         vhx(i) = vhx(i) + e2 * vx
         extot = extot + e2 * ex * rho(i) * fpi * r2(i) * r(i) * dx
         vxtot = vxtot + e2 * vx * rho(i) * fpi * r2(i) * r(i) * dx
      end do

      return
      end 

c--------------------------------------------------------------------
      subroutine slaterx(rho,ex,vx)
c--------------------------------------------------------------------
c Exchange energy per particle and potential in the Slater approx.
c J.C.Slater, Phys.Rev. 81, 385 (1951).
c The Xalpha scheme with alpha=2/3 (the exact DFT value) is used 
c W.Kohn and L.J.Sham Phys.Rev. 137,A1441 (1965).
c
      implicit none
      real*8 pi, third, alpha
      parameter (pi=3.14159265358979d0, third=1.d0/3.d0,
     +           alpha=2.d0/3.d0)
      real*8 rho, ex, vx, fac

      fac =-9.d0/8.d0*alpha*(3.d0/pi)**third

      ex = fac*rho**third
      vx = 4.d0/3.d0*ex

      return
      end
