program atom
      implicit double precision (a-h,o-z)
      dimension dval(2,561),nval(2,7,10),ddval(2,561),ndval(2,7,10)
      dimension dcore(2,561),vhv(561),fpu(561)
c
c free atom program  by m.p. summer 1985 (valence energy and core
c density are calculated)
c
c GGA implemented by T. Holmquist and U. Yxklinten.
c 
c files used:
c               6: terminal output
c              22: density and potential output
c              33: density and potential input
c              44: potential and density double prec. output
c              46: potential and density double prec. input
c              55: control input
c              66: printer output
c              77: density (total and core) output (J.H.)
c
      dimension r(561),veff(2,561),vc(561),vcold(561),vv(561)
      dimension gr(561),lm(2),rnocc(2,7,10),ux1(561),ux2(561)
      dimension u(561),roo(561),spl(561),nb(7)
      dimension db(2,561),v(561),vold(2,561)
      dimension snlo(561),wj(301),ekin(2)
      dimension ebound(2,7,10),nbound(2,7),de(2,7,10)
      dimension apu(561),bpu(561),cpu(561),dpu(561),nlp(2,7,10)
      dimension epu(561)
c===================================================================
c  gga
      dimension excgga(561),uxcup(561),uxcdn(561)
      logical doned,donela,donegr,donez,donec
c  gga
c===================================================================
      common/sc/gr,r,snlo,nbl
      common/mess/wj,dx,nblock,jblock
      common/pot/v
      frs(x)=(3.d0/(4.d0*pi*x))**(1.d0/3.d0)
      pi=4.d0*datan(1.d0)
      pi2=pi/2.d0
      a=(4.d0/9.d0/pi)**(1.d0/3.d0)
c
      open(22,file='adensout.dat',status='unknown')
      open(33,file='adensin.dat',status='unknown')
      open(44,file='apotout.dat',status='unknown')
      open(46,file='apotin.dat',status='unknown')
      open(55,file='atomctrl.dat',status='old')
      open(66,file='aprtout.dat',status='unknown')
      open(77,file='atomdens.dat',status='unknown')
      write(6,*) '=== Its ===== Energy (Ha) ======'
      call flush(6)
      read(55,*)z,zion
c     z     : atomic number
c     zion  : ionicity
c
      write(66,2)  z,zion
2     format(/' atom number:',f5.0,'  ionicity:',f5.0)
      sf=4.d0*pi
      read(55,*)jblock,nbl,c
c parameters of the hermann-skillmann mesh
      read(55,*)fback,thresh,qsc
C     fback   : feedback
c     thresh  : error allowed for the bound state eigenenergy
c               e.g. 0.00001
c     qsc     : parameter for the screened green's function (0.005)
      read(55,*)itmax,inopt,ipr,iout
c     itmax  : max. no. of iterations
c     iopt   : =0 starting potential generated; =1 starting potential
c               read in; =2 starting potential read in from unit 46.
c               (double precision format)
c     ipr, iout : print out parameters
c
c set up the herman-skillman mesh
c
      nblock=jblock
      mesh=nblock*nbl+1
      n=mesh
      mest=mesh
      dx=c*0.0025d0
      en0=-0.01d0
      i=1
      r(i)=0.d0
      deltax=dx
      do 251 j=1,nblock
      do 241 jk=1,nbl
      i=i+1
 241   r(i)=r(i-1)+deltax
      deltax=2.d0*deltax
  251  continue
c
      read(55,*)ne,iys,ixc
c     ne: no. bound states
c     iys: =1 for spin compensated; =2 spin polarized
c     ixc: =  1 b-h xc; =0 c-a xc
      if(ixc.eq.1)write(66,2149)
      if(ixc.eq.0)write(66,2148)
      if(ixc.eq.2)write(66,2147)
2149  format(/' von barth - hedin xc')
2148  format(/' ceperley - alder xc')
2147  format(/' perdew - wang xc')
      dxx=r(n)-r(n-1)
      write(66,5)jblock,nbl,c
      write(66,5849)dx,dxx,r(n)
  5   format(/' hermann-skillmann mesh',/,' blocks, nbl, c:',2(i5,','),
     j f10.7)
5849  format(' first interval, last interval, final r:',2(f10.6,','),
     j f10.6)
      write(66,2312) thresh,qsc,fback
2312  format(/' energy eigenvalue threshold, screening parameter, ',
     j 'feedback:',2(f12.6,','),f6.3)
      write(66,9)
      do 5273 i=1,7
5273  nb(i)=0
      lmax=0
9     format(//' initial bound state configuration',/,
     j '   nlm    energy    occup.  val.el. d-band el.')
      zv=0.d0
      do 5595 ispin=1,iys
      do 5595 ii=1,ne
      read(55,*)nnlz,eb0,rnoc,nv,ndv
c
c nnlz :  e.g. 100
c eb0  :   energy eigenvalue guess
c rnoc :   occupation number
c nv   :   =1 when orbital is a valence orbital; =0 for core orbitals
c ndv  :   =1 when orbital is a d-valence orbital; =0 for other orbitals
c
      nnn=nnlz/100
      lll=(nnlz-nnn*100)/10
      nnn=nnn-lll
      if(nnn.gt.nb(lll+1))nb(lll+1)=nnn
      if(lll.gt.lmax)lmax=lll
      rnocc(ispin,lll+1,nnn)=rnoc
        nval(ispin,lll+1,nnn)=nv
        ndval(ispin,lll+1,nnn)=ndv
        if(iys.eq.1)nval(2,lll+1,nnn)=nv
        if(iys.eq.1)ndval(2,lll+1,nnn)=ndv
          if(iys.eq.1)rnocc(2,lll+1,nnn)=rnoc
          if(iys.eq.1)rnoc=rnoc*2
      write(66,19)nnlz,eb0,rnoc,nv,ndv
      zv=zv+rnoc*(nv+ndv)
5595   ebound(ispin,lll+1,nnn)=eb0
       write(66,5120)zv
5120   format(/' number of valence electrons: ',f5.1)
 19   format(i6,e12.4,f7.2,i7,i10)
      nb1=nbl+1
      iter=0
      itr=itmax-iout
c
  16  format(7(f12.5,1x))
  17  format(/' charge and spin density profiles '/)
  18  format(1h0/' potentials vs. distance after ',i3,' iterations'/)
c
c tabulate screened green's function
      do 10 i=1,n
  10  gr(i)=exp(-qsc*r(i))
c
c integration weigths
      do 121 i=1,nb1
 121  wj(i)=(3.d0+(-1.d0)**i)/3.d0
      wj(1)=1.d0/3.d0
      wj(nb1)=1.d0/3.d0
c
c initial values of the potentials
      if(inopt.eq.1) goto 302
      if(inopt.eq.2) goto 306
c
      itt=0
c thomas-fermi potential
      do 30 i=2,n
      x=r(i)/0.88534135d0*z**(1.d0/3.d0)
      xx=sqrt(x)
      vc(i)=-z/r(i)/(1.+0.02747d0*xx+1.243d0*x-0.1486d0*x*xx
     j+0.2302d0*x*x+0.007298d0*x*x*xx+0.006944d0*x*x*x)
      do 30 ispin=1,2
 30   veff(ispin,i)=vc(i)
      goto 311
c     read in an old potential for the input
302   read(33,2645)itt
2645  format(i3)
1645  format(i3,' z,ion:',2f4.0,' jblock,nbl,c:',2i5,f10.7)
      do 303 i=1,n
 303  read(33,907) vc(i),veff(1,i),veff(2,i),du1,du2
      goto 311
 306  read(46,2645) itt
      do 304 i = 1, n
        read(46,*) du1, vc(i), vve, du2
        veff(1,i) = vve
        veff(2,i) = vve
 304  continue
 311  do 312 i=1,n
      vcold(i)=vc(i)
      do 312 ispin=1,2
 312  vold(ispin,i)=veff(ispin,i)
 5555 format(e15.8)
 7771 format(/' initial potentials'/)
 7773 format(5f13.5)
      if(iout.lt.-1)go to 1437
      write(66,7771)
      do 7772 i=2,n,ipr
 7772 write(66,7773) r(i),vc(i),veff(1,i),veff(2,i)
1437  eold=0.d0
      etot=-1.d0
      itec=0
c
c iteration
c
  100 iter=iter+1
c	write(6,*) 'NEW ITERATION'
c	call flush(6)
c convergency check
      if(abs(etot-eold).lt.5.d-5)itec=itec+1
      if(abs(etot-eold).lt.5.d-5.and.itec.eq.1)itmax=iter
      eold=etot
      if(iter-itmax)999,999,400
999   itt=itt+1
      write(66,6147)itt
      do 117 i=1,n
         dval(2,i)=0.d0
         ddval(2,i)=0.d0
         dcore(1,i)=0.d0
         dcore(2,i)=0.d0
         dval(1,i)=0.d0
 117     ddval(1,i)=0.d0
 6147    format(///' *********************   iteration:',
     1        i4,'  ***********************'/)
      do 917 ispin=1,iys
      do 116 i=1,n
      vv(i)=veff(ispin,i)
      v(i)=2*vv(i)
 116  db(ispin,i)=0.d0
      lm(ispin)=-1
      do 807 i=0,4
                ik=i
      if(i)817,817,818
817   do 819 jj=1,4
819   u(jj)=r(jj)-z*r(jj)**2
      go to 822
818   do 821 jj=1,4
821   u(jj)=r(jj)**(i+1)
822   call schrhs(vv,0.d0,ik,u)
      merkki=1
c determine the number of bound states
      ncross=0
      do 826 ii=2,n
      if(merkki)823,823,824
823   if(u(ii))826,826,825
824   if(u(ii))825,826,826
825   merkki=-merkki
      ncross=ncross+1
826   continue
      dlo=(u(n)-u(n-1))*u(n)
      if(dlo.ge.0)nbound(ispin,i+1)=ncross
      if(dlo.lt.0)nbound(ispin,i+1)=ncross+1
      if(nbound(ispin,i+1).eq.0)go to 827
      lm(ispin)=i
      write(66,678)i,nbound(ispin,i+1)
      nbound(ispin,i+1)=min0(nbound(ispin,i+1),nb(i+1))
678   format(' l = ',i4,',',i4,' bound states ')
807   continue
827   if(lm(ispin).lt.0.and.ispin.eq.2)go to 674
      if(lm(ispin).lt.0.and.ispin.eq.1)go to 917
      lm(ispin)=min0(lm(ispin),lmax)
      lmm=lm(ispin)
      do 918 i=0,lmm
      nii=min0(nbound(ispin,i+1),8)
                ik=i
      do 918 ii=1,nii
      if(rnocc(ispin,i+1,ii).lt.0.1d0)go to 918
      nn=i+ii
      en=2.d0*ebound(ispin,i+1,ii)
      if(en.ge.0.d0)en=en0
      en1=en
      dde=de(ispin,i+1,ii)
      if(iter.le.2)dde=0.d0
c     determine the bound state energy and eigenfunction
      call scheq(z,en,ik,nn,mest,mesh,c,thresh,iflag,npr,dde)
c scheq operates in rydberg units
      if(iter.ge.2)de(ispin,i+1,ii)=abs(en-en1)
      if(iflag.eq.1) goto 400
 163  ebound(ispin,i+1,ii)=en/2.d0
      nlp(ispin,i+1,ii)=npr
      do 9188 ji=1,n
        dval(ispin,ji)=dval(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/sf
     j*nval(ispin,i+1,ii)
        ddval(ispin,ji)=ddval(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/
     j sf*ndval(ispin,i+1,ii)
        dcore(ispin,ji)=dcore(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/
     j sf*(1-ndval(ispin,i+1,ii))*(1-nval(ispin,i+1,ii))
 9188  db(ispin,ji)=db(ispin,ji)+rnocc(ispin,i+1,ii)*snlo(ji)**2/sf
 918  continue
      do 919 ji=2,n
      dcore(ispin,ji)=dcore(ispin,ji)/r(ji)**2
 919  db(ispin,ji)=db(ispin,ji)/r(ji)**2
      db(ispin,1)=db(ispin,2)
      dcore(ispin,1)=dcore(ispin,2)
        do 5739 ji=2,n
        ddval(ispin,ji)=ddval(ispin,ji)/r(ji)**2
5739    dval(ispin,ji)=dval(ispin,ji)/r(ji)**2
        dval(ispin,1)=dval(ispin,2)
        ddval(ispin,1)=ddval(ispin,2)
 917  continue
      if(iys.eq.2)go to 923
      do 921 i=1,n
      dval(2,i)=dval(1,i)
      ddval(2,i)=ddval(1,i)
      dcore(2,i)=dcore(1,i)
921   db(2,i)=db(1,i)
      lmm=lm(1)
          lm(2)=lm(1)
      do 922 i=0,lmm
      nii=nbound(1,i+1)
          nbound(2,i+1)=nii
      do 922 ii=1,nii
922   ebound(2,i+1,ii)=ebound(1,i+1,ii)
923   continue
 8    format(' ',i6,' energy: ',f12.5,' iteration loops: ',i4)
c
c total bound state energy
c
c	write(6,*) 'Total energy starts'
c	call flush(6)
  661 eb=0.d0
      ebv=0.d0
      write(66,7496)
7496  format(//' bound states'/)
      do 673 ispin=1,2
      lmm=lm(ispin)
      do 673 i=0,lmm
      nii=nbound(ispin,i+1)
      do 673 j=1,nii
      nn=i+j
      nnlz=100*nn+10*i
      eb=eb+ebound(ispin,i+1,j)*rnocc(ispin,i+1,j)
      ebv=ebv+ebound(ispin,i+1,j)*rnocc(ispin,i+1,j)*
     j (nval(ispin,i+1,j)+ndval(ispin,i+1,j))
          if(rnocc(ispin,i+1,j).lt.0.1d0) go to 673
      write(66,8) nnlz,ebound(ispin,i+1,j),nlp(ispin,i+1,j)
 673  continue
      go to 679
674   write(66,676)
676   format('0no bound states')
679   continue
c
c total charge and spin densities
c computing energy integrals
c
506   continue
c=========================================
c  gga 
c  set the "done" controle variables to false in the beginning of
c  each iteration

      if (ixc.eq.2) then
         doned = .false.
         donela = .false.
         donegr = .false.
         donez = .false.
         donec = .false.
      endif
c  gga
c========================================
      do 660 i=2,n
      roo(i)=db(1,i)+db(2,i)
      rr=roo(i)
c When db = (0,0) then the spin-polarization, spl = 0, and not 0/0.
      if (rr.gt.0.0000000001) then
         spl(i)=(db(1,i)-db(2,i))/rr
      else
         spl(i)=0.0
      endif
      x1=r(i)**2
      apu(i)=x1*veff(1,i)*db(1,i)*sf
      bpu(i)=x1*veff(2,i)*db(2,i)*sf
      if (ixc.eq.2) then
c================================================================
c  gga
c  calculate the gga exchange-correlation energy

         call ggaexc(r,db,i,1,n,doned,excgga(i))
         cpu(i) = x1*rr*excgga(i)*sf
c  gga
c================================================================
      else
         cpu(i)=x1*rr*exc(rr,spl(i),ixc)*sf
      endif
      dpu(i)=x1*spl(i)*rr*sf
 660  continue
 26   format(' induced moment: ',f10.5)
      call simpsh(apu,v1)
      call simpsh(bpu,v2)
      call simpsh(cpu,eexc)
      call simpsh(dpu,smom)
c
c compute coulomb energy
c
      cz=z
      do 683 i=1,n
        bpu(i)=roo(i)*r(i)**2*sf
        dpu(i)=(dval(1,i)+dval(2,i))*r(i)**2*sf
        apu(i)=(ddval(1,i)+ddval(2,i))*r(i)**2*sf
	cpu(i)=roo(i)*(r(i)**2*vc(i)-z*r(i))*sf/2.d0
683   continue
      call simpsh(dpu,sum1)
      call simpsh(apu,sum3)
      call simpsh(bpu,sum2)
      call simpsh(cpu,ec)
  24  format(' kin: ',2(e14.6,1x),' coul: ',e14.6,' exc: ',e14.6)
  25  format(/' total energy: ',e15.7)
      write(66,192) sum2,sum1,sum3
      write(66,26) smom
      ekin(1)=-v1
      ekin(2)=-v2
      etot=ekin(1)+ekin(2)+eb+ec+eexc
      write(66,7453)
7453  format(/' energy terms')
      write(66,622)eb
 622  format(' energy eigenvalue sum: ',e14.6)
      write(66,24) ekin(1),ekin(2),ec,eexc
      write(66,25) etot
c      write(6,2573) itt,sum2,etot
c      call flush(6)
c2573  format(' iter:',i4,' total ch:',f6.2,' total en:',f15.7)
      write(6,2573) itt,etot
      call flush(6)
2573  format('    ',i4,'   ',f15.5)

      if(iter.ne.1 .and. iter.lt.itr) goto 167
      if(iout.lt.0)go to 167
      write(66,17)
  192  format(/' integr charges: total  ',f10.7,' valence: ',2f10.7)
      do 166 i=1,n,ipr
      r22=r(i)*r(i)
      db1=db(1,i)
      db2=db(2,i)
      dr=db1+db2
      dcc=dcore(1,i)+dcore(2,i)
      dvv=dval(1,i)+dval(2,i)+ddval(1,i)+ddval(2,i)
      ddb=db1-db2
c166  write(66,1624) r(i),db1,db2,dval(i),ddval(i),dr,ddb
 166  write(66,1624) r(i),db1,db2,dr,dcc,dvv
c
c
c compute coulomb potential
c
 167  do 55 i=2,n
  55  vc(i)=r(i)*vc(i)+zion
c	write(6,*) 'Coulomb potential'
c	call flush(6)
      zi=z-zion
      vc(1)=-zi
      call scrhs(roo,vc,zi,qsc,n)
      do 50 i=2,n
  50  vc(i)=(vc(i)-zion)/r(i)
c
c set up the total potential
c
      do 60 i=2,n
      dr=roo(i)
      vc(i)=(1.d0-fback)*vcold(i)+fback*vc(i)
      vcold(i)=vc(i)
      if (ixc.eq.2) then
c================================================================
c  gga
c  calculate the gga exchange-correlation potential

         call ggauxc(r,db,spl,i,1,n,donela,donegr,donez,
     j        uxcup(i),uxcdn(i))
         ux1(i) = uxcup(i)
         ux2(i) = uxcdn(i)
c  gga
c================================================================
      else   
         ux1(i)=uxc(dr,spl(i),1,ixc)
         ux2(i)=uxc(dr,spl(i),2,ixc)
      endif
 6848 veff(1,i)=vc(i)+ux1(i)
      veff(2,i)=vc(i)+ux2(i)
      do 60 ispin=1,2
      veff(ispin,i)=(1.d0-fback)*vold(ispin,i)+fback*veff(ispin,i)
  60  vold(ispin,i)=veff(ispin,i)
c  calculate the valence energy
      do 5100 i=2,mesh
      do 5110 j=1,mesh
      apu(j)=0.d0
      if(i.le.j)apu(j)=(dval(1,j)+dval(2,j)+ddval(1,j)+ddval(2,j))*
     j (r(j)**2/r(i)-r(j))
5110  continue
      call simpsh(apu,vp)
5100  vhv(i)=-vp*4.d0*pi+zv/r(i)
      do 5130 i=1,mesh
      x1=r(i)**2
      rr=db(1,i)+db(2,i)
      spt=0.d0
      if(rr.gt.1.d-10)spt=(db(1,i)-db(2,i))/rr
      rrv=dval(1,i)+dval(2,i)+ddval(1,i)+ddval(2,i)
      spv=0.d0
      if(rrv.gt.1.d-10)spv=(dval(1,i)-dval(2,i)+ddval(1,i)-ddval(2,i))
     j /rrv
      rrc=dcore(1,i)+dcore(2,i)
      spc=0.d0
      if(rrc.gt.1.d-10)spc=(dcore(1,i)-dcore(2,i))/rrc
      apu(i)=-0.5d0*x1*sf*rrv*vhv(i)
      if (ixc.eq.2) then
c=============================================================
c  gga
c  calculate the gga exchange-correlation energy for core electrons

         if (rrc.lt.1.d-100) then
            excgac = 0.d0
         else
            call ggaexc(r,dcore,i,1,n,donec,excgac)
         endif
         bpu(i)=x1*rr*excgga(i)*sf
         cpu(i)=-x1*rrc*excgac*sf
         fpu(i)=(dcore(1,i)+dcore(2,i))*excgga(i)*x1*sf
         dpu(i)=(-(dval(1,i)+ddval(1,i))*uxcup(i)
     j        -(dval(2,i)+ddval(2,i))*uxcdn(i))*sf*x1
c  gga
c==================================================================
      else
         bpu(i)=x1*rr*exc(rr,spt,ixc)*sf
         cpu(i)=-x1*rrc*exc(rrc,spc,ixc)*sf
         fpu(i)=(dcore(1,i)+dcore(2,i))*exc(rr,spt,ixc)*x1*sf
         dpu(i)=(-(dval(1,i)+ddval(1,i))*uxc(rr,spt,1,ixc)
     j        -(dval(2,i)+ddval(2,i))*uxc(rr,spt,2,ixc))*sf*x1
      endif
      epu(i)=apu(i)+bpu(i)+cpu(i)+dpu(i)+fpu(i)
5130  continue
c      call simpsh(apu,ec)
c      call simpsh(bpu,eext)
c      call simpsh(fpu,eexx)
c      call simpsh(cpu,eexc)
c      call simpsh(dpu,eexv)
c      etv=ebv+ec+eext+eexc+eexv
      call simpsh(epu,etv)
      etv = etv + ebv
      write(66,5140)etv
5140  format(' valence energy:',e15.7)
      if(iter.ne.1 .and. iter.lt.itr)  go to 100
      if(iout.lt.0)go to 100
      write(66,18) itt
      do 75 i=2,n,ipr
      dr=roo(i)
 75   write(66,1624) r(i),vc(i),ux1(i),ux2(i),veff(1,i),veff(2,i)
1624  format(1x,f12.5,6e14.5)
      goto 100
  400 continue
c
c
c
        do 1140 i=2,mesh
        do 1130 j=1,mesh
        apu(j)=0.d0
        if(i.le.j)apu(j)=(db(1,j)+db(2,j))*(r(j)**2/r(i)-r(j))
1130    continue
        call simpsh(apu,vp)
1140    vc(i)=-vp*sf
      write(22,1645)itt,z,zion,jblock,nbl,c
      write(44,1645)itt,z,zion,jblock,nbl,c
      do 900 i=1,n
          dd=db(1,i)+db(2,i)
          dc=dcore(1,i)+dcore(2,i)
	  dddc = dd + dc
          dv1=dval(1,i)+dval(2,i)
          ddv1=ddval(1,i)+ddval(2,i)
          vve = 0.5d00*(veff(1,i) + veff(2,i))
          write(44,*) r(i), vc(i), vve, dddc
c         The line below writes total and core densities to atomdens.dat (J.H.)
          write(77,906) r(i), 4*pi*dd*r(i)*r(i), 4*pi*dc*r(i)*r(i)
          write(22,907)vc(i),veff(1,i),veff(2,i),dddc,r(i)
 900  continue
 907  format(6e15.7)
c900  write(22,906)vc(i),veff(1,i),veff(2,i),dd,dc
c900  write(22,906)vc(i),veff(1,i),veff(2,i),dd,dv1,ddv1
 906  format(6e20.12)
      write(6,*) '================================'
      call flush(6)
c      ieeer=ieee_flags('clear','exception','all',ieeeout)
      close(6)
      close(22)
      close(33)
      close(55)
      close(66)
      close(77)
      stop
      end