MODULE AEatom
  USE GlobalMath
  USE excor
  USE radialsch
  USE calcpotential
  USE anderson_realmix
  USE atomdata

  IMPLICIT NONE
  REAL(8), PARAMETER, PRIVATE ::small0=1.d-15,rimix=0.5d0,worst=1.d-5
  REAL(8), PARAMETER, PRIVATE ::range0=50.d0
  INTEGER, PARAMETER, PRIVATE :: niter=1000,mxgrid=20000,mxloop=200
  REAL(8), PRIVATE :: conv1=1.d13,conv2=2.d13,conv3=3.d13,conv4=4.d13

CONTAINS

  SUBROUTINE SCFatom(AEGrid,AEPot,AEOrbit,AESCF)
     TYPE (GridInfo), INTENT(INOUT) :: AEGrid
     TYPE (PotentialInfo), INTENT(INOUT) :: AEPot
     TYPE (OrbitInfo), INTENT(INOUT) :: AEOrbit
     TYPE (SCFInfo), INTENT(INOUT) :: AESCF

    !  program to calculate the self-consistent density functional
    !    atom ground state fo atom with atomic number nz
    !  iscr1,iscr2 are scratch units
    !    for self-consistent potential rv

    TYPE (Anderson_Context) , POINTER :: AErho

    INTEGER :: nps,npp,npd,npf,npg,norbit,nz,n
    REAL(8) :: electrons,h
    REAL(8) :: rmix,xocc,qf,small,range,zeff,delta
    REAL(8) :: qcal, rescale,cnvrg,emin,ecoul,eexc,etxc,eone,etot
    INTEGER :: icount,i,j,k,it,start,np,ierr,nroot
    INTEGER :: is,ip,id,jf,ig,io,l,nfix,ir,nzeff,loop,jierr
    REAL(8), ALLOCATABLE :: denout(:)
    INTEGER, ALLOCATABLE :: nl(:,:)

    WRITE(6,*) 'enter atomic symbol and atomic number'
    READ(5,*) AEPot%sym,nz
    WRITE(6,*) 'exchange-correlation type -- LDA-PW(default) or GGA-PBE '
    READ(5,*) exctype
    AEPot%nz=nz
    WRITE(6,*) 'Calculation for atomic number = ',nz
    WRITE(6,*) 'enter maximum principal quantum numbers for s,p,d,f,g'
    READ(5,*) nps,npp,npd,npf,npg
    IF(nps<0)nps=0
    IF(npp<0)npp=0
    IF(npd<0)npd=0
    IF(npf<0)npf=0
    IF(npg<0)npg=0
    WRITE(6,'(5i4)') nps,npp,npd,npf,npg
    i=MAX(nps,npp,npd,npf,npg)
    j=nps+npp+npd+npf+npg
    ALLOCATE(nl(i,j),AEOrbit%np(j),AEOrbit%l(j),&
          AEOrbit%eig(j),AEOrbit%occ(j),stat=k)
    IF (k/=0) THEN
       WRITE(6,*) 'Error in allocation of nl,np...',k
       STOP
    ENDIF
     nl=0;AEOrbit%np=0;AEOrbit%l=0;AEOrbit%eig=0;AEOrbit%occ=0
    icount=0
    IF (nps.GT.0) THEN
       DO is=1,nps
          icount=icount+1
          nl(is,1)=icount
          AEOrbit%occ(icount)=2.d0
          AEOrbit%np(icount)=is
          AEOrbit%l(icount)=0
       ENDDO
    ENDIF
    IF (npp.GT.1) THEN
       DO ip=2,npp
          icount=icount+1
          nl(ip,2)=icount
          AEOrbit%occ(icount)=6.d0
          AEOrbit%np(icount)=ip
          AEOrbit%l(icount)=1
       ENDDO
    ENDIF
    IF (npd.GT.2) THEN
       DO id=3,npd
          icount=icount+1
          nl(id,3)=icount
          AEOrbit%occ(icount)=10.d0
          AEOrbit%np(icount)=id
          AEOrbit%l(icount)=2
       ENDDO
    ENDIF
    IF (npf.GT.3) THEN
       DO jf=4,npf
          icount=icount+1
          nl(jf,4)=icount
          AEOrbit%occ(icount)=14.d0
          AEOrbit%np(icount)=jf
          AEOrbit%l(icount)=3
       ENDDO
    ENDIF
    IF(npg.GT.4) THEN
       DO ig=5,npg
          icount=icount+1
          nl(ig,5)=icount
          AEOrbit%occ(icount)=18.d0
          AEOrbit%np(icount)=ig
          AEOrbit%l(icount)=4
       ENDDO
    ENDIF
    norbit=icount
    AEOrbit%nps=nps
    AEOrbit%npp=npp
    AEOrbit%npd=npd
    AEOrbit%npf=npf
    AEOrbit%npg=npg
    AEOrbit%norbit=norbit
    WRITE(6,*) norbit, ' orbitals will be calculated'
    !
    WRITE(6,*)' Below are list the default occupations '
    WRITE(6,"(' n  l     occupancy')")
    DO io=1,norbit
       WRITE(6,'(i2,1x,i2,4x,1pe15.7)') &
          AEOrbit%np(io),AEOrbit%l(io),AEOrbit%occ(io)
    ENDDO
    !
    WRITE(6,*)' enter np l occ for all occupations for all revisions'
    WRITE(6,*)'  enter 0 0 0. to end'

    DO
       READ(5,*) ip,l,xocc
       IF (ip.LE.0) EXIT
       nfix=nl(ip,l+1)
       IF (nfix.LE.0.OR.nfix.GT.norbit) THEN
          WRITE(6,*) 'error in occupations -- ip,l,xocc',                  &
               ip,l,xocc,nfix,norbit
          STOP
       ENDIF
       AEOrbit%occ(nfix)=xocc
    ENDDO

    !
    WRITE(6,*) ' Corrected occupations are: '
    WRITE(6,"(' n  l     occupancy')")
    electrons=0.d0
    DO io=1,norbit
       WRITE(6,'(i2,1x,i2,4x,1pe15.7)')&
           AEOrbit%np(io),AEOrbit%l(io),AEOrbit%occ(io)
       electrons=electrons+AEOrbit%occ(io)
    ENDDO
    AEPot%q=electrons
    qf=nz-electrons
    WRITE(6,*)
    WRITE(6,*) 'nuclear charge    = ' , nz
    WRITE(6,*) 'electronic charge = ', electrons
    WRITE(6,*) 'net charge        = ', qf
    !
    n=mxgrid

    ALLOCATE(AEGrid%r(n),AEPot%den(n),AEPot%rv(n),&
           AEOrbit%wfn(n,j),denout(n),stat=k)
    IF (k/=0) THEN
       WRITE(6,*) 'Error in allocation of r,den,...',k
       STOP
    ENDIF
    small=small0
    rmix=rimix
    range=range0

    h=range/n
    WRITE(6,*) 'n range h = ', n,range,h
    AEGrid%n=n
    AEGrid%h=h
    DO ir=1,n
       AEGrid%r(ir)=ir*h
    ENDDO
    IF (small.LT.0.d0.OR.small.GT.1.d0) small=small0

    !
    !  calculate initial charge density from hydrogen-like functions
    !    also initial energies
    !
    AEPot%den(1:n)=0.d0
    zeff=nz
    DO io=1,norbit
       np=AEOrbit%np(io)
       l=AEOrbit%l(io)
       xocc=AEOrbit%occ(io)
       zeff=zeff-0.5d0*xocc
       zeff=MAX(zeff,1.d0)
       nzeff=zeff+0.1d0
       AEOrbit%eig(io)=-(zeff/(np))**2
       WRITE(6,*) io,np,l,xocc,AEOrbit%eig(io)
       DO ir=1,n
          AEOrbit%wfn(ir,io)=hwfn(nzeff,np,l,AEGrid%r(ir))
          AEPot%den(ir)=AEPot%den(ir)+xocc*(AEOrbit%wfn(ir,io)**2)
       ENDDO
       zeff=zeff-0.5d0*xocc
    ENDDO
    !
    !  check charge
    !
    qcal=overint(n,h,AEPOT%den)
    qf=qcal
    WRITE(6,*) 'qcal electrons = ',qcal, electrons
    !  rescale density
    rescale=electrons/qcal
    AEPot%den(1:n)=AEPot%den(1:n)*rescale 
    !
    !  choose form of exchange-correlation potential
    CALL initexch
    !
    cnvrg=worst
    WRITE(6,*) 'Density convergence parameter set to at least',cnvrg
    !
    !    start iteration loop
    !
    DO   loop=1,mxloop
       !  calculate  potential * r == rv for given density
       !
       CALL potential(AEGrid,AEPot,ecoul,etxc,eexc)
       !
       !  solve for bound states of Schroedinger equation
       !
       icount=0
       qf=nz-electrons
       jierr=0
       it=0
       !  s states :
       IF (nps.GT.0) THEN
          it=it+1
          emin=-nz*nz
          l=0
          nroot=nps
          start=1
          CALL boundsch(AEGrid,AEPot,AEOrbit,l,start,nroot,emin,ierr)
       ENDIF
       !  p states :
       IF (npp.GT.1) THEN
          it=it+1
          emin=-nz*nz/4.d0
          l=1
          nroot=npp-1
          start=start+nps
          CALL boundsch(AEGrid,AEPot,AEOrbit,l,start,nroot,emin,ierr)
       ENDIF
       !  d states :
       IF (npd.GT.2) THEN
          it=it+1
          emin=-nz*nz/9.d0
          l=2
          nroot=npd-2
          start=start+npp-1
          CALL boundsch(AEGrid,AEPot,AEOrbit,l,start,nroot,emin,ierr)
       ENDIF
       !  f states :
       IF (npf.GT.3) THEN
          it=it+1
          emin=-nz*nz/16.d0
          l=3
          nroot=npf-3
          start=start+npd-2
          CALL boundsch(AEGrid,AEPot,AEOrbit,l,start,nroot,emin,ierr)
       ENDIF
       !  g states :
       IF (npg.GT.4) THEN
          it=it+1
          emin=-nz*nz/25.d0
          l=4
          nroot=npg-4
          start=start+npf-3
          CALL boundsch(AEGrid,AEPot,AEOrbit,l,start,nroot,emin,ierr)
       ENDIF
       !
       !  calculate new density
       !
       denout(1:n)=0.d0 
       WRITE(6,*) '  results for loop = ',loop
       WRITE(6,*) ' n  l     occupancy       energy'
       DO io=1,norbit
          WRITE(6,'(i2,1x,i2,4x,1p2e15.7)') &
               AEOrbit%np(io),AEOrbit%l(io),&
               AEOrbit%occ(io),AEOrbit%eig(io)
          IF (AEOrbit%occ(io).GT.small) THEN
             DO i=1,n
                denout(i)=denout(i)+AEOrbit%occ(io)*(AEOrbit%wfn(i,io)**2)
             ENDDO
          ENDIF
       ENDDO
       qcal=overint(n,h,denout)
       WRITE(6,*) 'qcal electrons = ',qcal, electrons
       !  rescale density
       rescale=electrons/qcal
       denout(1:n)=denout(1:n)*rescale
       !
       denout=denout-AEPot%den
       delta=SUM(ABS(denout))
       CALL shift4(conv1,conv2,conv3,conv4,delta)
       WRITE(6,*) 'density iter',loop,delta
       IF (loop.EQ.1) CALL InitAnderson(AErho,6,5,n,rmix,1.d5) 
       CALL Anderson_Mix(AErho,AEPot%den(1:n),denout(1:n))
       !  correct and recale density
       DO i=1,n
          IF (AEPot%den(i).LT.0.d0) AEPot%den(i)=0.d0
       ENDDO
       qcal=overint(n,h,AEPot%den)
       WRITE(6,*) 'qcal electrons = ',qcal, electrons
       !  rescale density
       rescale=electrons/qcal
       AEPot%den(1:n)=AEPot%den(1:n)*rescale 
       !
       WRITE(6,*) '  results for loop ,delta = ',loop,delta
       WRITE(6,*) ' n  l     occupancy       energy'
       eone=0.d0
       DO io=1,norbit
          WRITE(6,'(i2,1x,i2,4x,1p2e15.7)') &
               AEOrbit%np(io),AEOrbit%l(io),&
               AEOrbit%occ(io),AEOrbit%eig(io)
          eone=eone+AEOrbit%occ(io)*AEOrbit%eig(io)
       ENDDO
       WRITE(6,*)
       WRITE(6,*) ' Total energies'
       WRITE(6,*) '    One-electron contribution:  ',eone
       WRITE(6,*) '    Coulomb contribution     :  ',ecoul
       WRITE(6,*) '    Exch-correl contribution :  ',eexc
       etot=eone-ecoul+etxc
       WRITE(6,*) '    Total                    :  ',etot
       !
       IF (loop>=4) THEN
          IF (.NOT.(conv4.LE.conv3.AND.conv3.LE.conv2 &
               .AND.conv2.LE.conv1).AND.conv4.LE.cnvrg) THEN
             !
             !  converged result
             !
             WRITE(6,*) '  dfatom converged in',loop,' iterations'
             AESCF%iter=loop
             WRITE(6,*) '     for nz = ',nz
             WRITE(6,*) '    delta(density)  = ', delta
             AESCF%delta=delta
             WRITE(6,*) '  results for loop = ',loop
             WRITE(6,*) ' n  l     occupancy       energy'
             eone=0.d0
             DO io=1,norbit
                WRITE(6,'(i2,1x,i2,4x,1p2e15.7)') &
                     AEOrbit%np(io),AEOrbit%l(io),&
                     AEOrbit%occ(io),AEOrbit%eig(io)
                eone=eone+AEOrbit%occ(io)*AEOrbit%eig(io)
             ENDDO
             WRITE(6,*)
             WRITE(6,*) ' Total energies'
             WRITE(6,*) '    One-electron contribution:  ',eone
             AESCF%eone=eone
             WRITE(6,*) '    Coulomb contribution     :  ',ecoul
             AESCF%ecoul=ecoul
             WRITE(6,*) '    Exch-correl contribution :  ',eexc
             AESCF%eexc=eexc
             etot=eone-ecoul+etxc
             AESCF%etot=etot
             WRITE(6,*) '    Total                    :  ',etot
             CALL FreeAnderson(AErho)
             DEALLOCATE(denout,nl)
             RETURN
          ENDIF
       ENDIF
    ENDDO   ! mxloop

    WRITE(6,*)'calculation terminating without density convergence'
    WRITE(6,*) 'delta, cnvrg =', delta,cnvrg
    WRITE(6,*) 'loop,mxloop = ',loop,mxloop
    CALL FreeAnderson(AErho)
    DEALLOCATE(denout,nl)
    STOP
  END SUBROUTINE SCFatom

  FUNCTION factorial(n)
    REAL(8) :: factorial
    INTEGER, INTENT(IN) :: n
    INTEGER :: i
    factorial=1
    IF (n.LT.2) RETURN
    DO i=2,n
       factorial=factorial*i
    ENDDO
  END FUNCTION factorial
  FUNCTION hwfn(nz,np,l,r)
    ! function to calculate the radial H wfn for nuclear charge nz
    !                                            principal qn   np
    !                                            orbital qn     l
    !   r*(radial H wfn) is returned
    REAL(8) :: hwfn
    REAL(8), INTENT(IN) :: r
    INTEGER, INTENT(IN) :: nz,np,l
    INTEGER :: node,k
    REAL(8) :: scale,rho,pref,term,sum
    node=np-l-1
    scale=2.d0*nz/np
    rho=scale*r
    pref=scale*SQRT(scale*factorial(np+l)/(2*np*factorial(node)))
    term=(rho**l)/factorial(2*l+1)
    sum=term
    IF (node.GT.0) THEN
       DO k=1,node
          term=-term*(node-k+1)*rho/(k*(2*l+1+k))
          sum=sum+term
       ENDDO
    ENDIF
    hwfn=r*pref*EXP(-0.5d0*rho)*sum
    !     write(6,*) 'r,hwfn=',r,hwfn
  END FUNCTION hwfn

END MODULE AEatom