onlineJ.html

c     Hedin-Lundqvist version 2001 April, EHCC, Kyoto University
c     -----non-relativistic LSD atomic structure calculation---------
c          writen by Masahiro Yamamoto,
c                 Institute of Atomic Enenrgy, Kyoto University
c                   March 1994
c     ============== INPUT ===============
c       file:atomlsd.dat
c             z : atomic number
c             no: # of orbitals considered
c                 spin up and down are different orbitals
c             npq(no): pricipal quantum number
c             elu(no): angular momentum quantum number
c             nspin(no): spin quantum number up=1 down=2
c             wei(no) : # of electrons
c             e(no)   : orbital energy in au
c             rv(2,651) : spin-dependent potential r*v(r)_sigma
c             iskip: if iskip=1,eigenvalue and eigenfunction are
c               calculated for wei=0.0
c     ============ OUTPUT ===================
c     FILE='RSCFEN.DAT':eigenvalue
c     FILE='RSCFPO.DAT':potential for up spin and down spin
c     FILE='ROSUM.DAT':electron density,spin up and down density
c     FILE='KINETIC.DAT':kinetic energy
c     file='wavef.dat':wavefunctions
      implicit real*8 (a-g,o-z)
      dimension npq(40),elu(40),nspin(40),wei(40),e(40)
      dimension rv(2,651),r(651),x(651),f(651),g(651),rhs(651)
      dimension xx(4),y(4),s(4),cof(4),a(5),b(4),p1(2,651),p2(2,651)
      dimension p3(2,651),fo(651),go(651),fw(40,651),rosum(2,651)
      dimension ro(651),rvsig(2,651),ek(40)
c------------------- data input ------------------
      open (1,file='atomlsd.dat')
      read (1,*) z
      read (1,*) no,iskip
      write (6,*) 'z=',z,'no=',no
      do i =1 , no
      read (1,*) npq(i),elu(i),nspin(i),wei(i),e(i)
      write (6,*) npq(i),elu(i),nspin(i),wei(i),e(i)
      enddo
      read (1,20) (rv(1,i), i=1, 651)
      read (1,20) (rv(2,i), i=1, 651)
   20 format (5(1pd15.7))
      close (1)
      write (6,20) (rv(1,i), i=1, 20)
      write (6,20) (rv(2,i), i=1, 20)
c     -----coordinate----log. mesh --------------------------
      dx=1.0d0/50.0d0
      do 40 i=1, 651
       x(i)=-9.0d0+dble(i-1)*dx
       r(i)=dexp(x(i))
   40 continue
c     ---------------potential-----------------------
      do i=1 , 651
       do j=1,2
       p1(j,i)=rv(j,i)
       p2(j,i)=rv(j,i)
       p3(j,i)=rv(j,i)
       enddo
      enddo
c     ------scf loop 1000-------------------------------
      do 1000 iscf=1, 1000
c     -------------- nlspin loop 2000----------
       do 2000 inls=1, no
        if ( iskip .eq. 1) go to 52
        if (wei(inls) .le. 1.0d-4) go to 2000
   52     iflag=0
c       --------- classical turning point ---------------
2010   iflag=iflag+1
        do 55 i=1, 650
         f(i)=0.0
         g(i)=0.0
         fo(i)=0.0
         go(i)=0.0
   55   continue
        do 50 i= 1, 651
         if ((rv(nspin(inls),i)/r(i)-e(inls)) .ge. 0.0) go to 60
   50   continue
   60   im=i
        if (mod(im,2) .eq. 0) im=im+1
c       -------- outermost point ---------
        do 70 i=im, 651
         if ((rv(nspin(inls),i)/r(i)-e(inls))*r(i)*r(i) .ge. 75.0)
     &    go to 80
   70   continue
   80   it=i
        if (mod(it,2) .eq. 0) it=it+1
        if (it .ge. 651) it=651
        write (6,*) 'classical t. p.', im,'outermost point',it
c     --------------------- radial sch. eq. r.h.s. ---------
      do 110 j=1, 651
         rhs(j)=2.0d0*rv(nspin(inls),j)-2.0d0*e(inls)*r(j)
         if (nint(elu(inls)) .eq. 0) go to 110
         rhs(j)=rhs(j)+(elu(inls)*(elu(inls)+1.0d0))/r(j)
110   continue
c      ------------- r:small region / power expansion approxi. -----
c      --- potential rv is expanded by sum[i=0to3]b(i)*r**i near r=0
       do 112 i=1, 4
        y(i)=rv(nspin(inls),i)
        xx(i)=r(i)
112 continue
      n=4
      do 11 i=1, n
       s(i)=0.
       cof(i)=0.
11 continue
      s(n)=-xx(1)
      do 13 i=2, n
       do 12 j=n+1-i, n-1
        s(j)=s(j)-xx(i)*s(j+1)
12   continue
       s(n)=s(n)-xx(i)
13 continue
      do 16 j=1,n
       phi=n
       do 14 k=n-1, 1, -1
        phi=k*s(k+1)+xx(j)*phi
14   continue
       fff=y(j)/phi
       bb=1.
       do 15 k=n,1,-1
        cof(k)=cof(k)+bb*fff
        bb=s(k)+xx(j)*bb
15   continue
16 continue
c      write (6,*) (i,cof(i), i=1, n)
      do 113 i=1, 4
       b(i)=cof(i)
113 continue
c     --near r=0 wave function f is expanded by sum[i=0to4]a(i)*r**i
      eluu=elu(inls)
      ee=e(inls)
      a(1)=1.0d0
      a(2)=b(1)*a(1)/(eluu+1.0d0)
      a(3)=(b(1)*a(2)+(b(2)-ee)*a(1))/(2.0d0*eluu+3.0d0)
      a(4)=(b(1)*a(3)+(b(2)-ee)*a(2)+b(3)*a(1))/(3.0d0*eluu+6.0d0)
      a(5)=(b(1)*a(4)+(b(2)-ee)*a(3)+b(3)*a(2)+b(4)*a(1))
     1     /(4.0d0*eluu+10.0d0)
      do 116 j=1, 4
       f(j)=0.0
       g(j)=0.0
       do 114 i=1, 5
        iram=i-1
        f(j)=f(j)+a(i)*r(j)**(nint(eluu)+1+iram)
        g(j)=g(j)+a(i)*(eluu+1.0+dble(iram))*r(j)**(nint(eluu)+iram)
114   continue
116 continue
C      ***** OUTWARD ADAMS-MOULTON METHOD R(4)->R(IM-) ************
       I=4
       FF1=DEXP(X(I))*G(I)
       FF2=DEXP(X(I-1))*G(I-1)
       FF3=DEXP(X(I-2))*G(I-2)
       FF4=DEXP(X(I-3))*G(I-3)
       GG1=rhs(i)*F(I)
       GG2=rhs(i-1)*F(I-1)
       GG3=rhs(i-2)*F(I-2)
       GG4=rhs(i-3)*F(I-3)
C      ----------- OUTWARD LOOP 1001 ----------
       DO 1001 I=4, im-1
C       ------- PREDICT -------
        FPR=F(I)+dx/24.0d0*(55.0d0*FF1-59.0d0*FF2+37.0d0*FF3-9.0d0*FF4)
        GPR=G(I)+dx/24.0d0*(55.0d0*GG1-59.0d0*GG2+37.0d0*GG3-9.0d0*GG4)
C       ------- CORRECT -------
        FF0=DEXP(X(I+1))*GPR
        GG0=rhs(i+1)*FPR
        FCR=F(I)+dx/24.0d0*(9.0d0*FF0+19.0d0*FF1-5.0d0*FF2+FF3)
        GCR=G(I)+dx/24.0d0*(9.0d0*GG0+19.0d0*GG1-5.0d0*GG2+GG3)
        DO 111 J=1, 2
         FF0=DEXP(X(I+1))*GCR
         GG0=rhs(I+1)*FCR
         FCR=F(I)+dx/24.0d0*(9.0d0*FF0+19.0d0*FF1-5.0d0*FF2+FF3)
         GCR=G(I)+dx/24.0d0*(9.0d0*GG0+19.0d0*GG1-5.0d0*GG2+GG3)
111   CONTINUE
        F(I+1)=FCR
        G(I+1)=GCR
C       ------------ PARAMETER EXCHANGE FOR THE NEXT LOOP -------
        FF4=FF3
        FF3=FF2
        FF2=FF1
        FF1=DEXP(X(I+1))*G(I+1)
        GG4=GG3
        GG3=GG2
        GG2=GG1
        GG1=rhs(i+1)*F(I+1)
1001 CONTINUE
c     ----check of the node ------------------------------
      iinod=0
      do 150 i=1, im-1
       if (f(i)*f(i+1) .le. 0.0) iinod=iinod+1
150 continue
      inode=npq(inls)-nint(elu(inls))-1
      if (inode-iinod) 155,170,160
155 e(inls)=1.25*e(inls)
      write (6,*) '# of node',inode,'# of node in this cal',iinod
      go to 2010
160 e(inls)=0.75*e(inls)
      write (6,*) '# of node',inode,'# of node in this cal',iinod
      go to 2010
c -------- inward adams-moulton method r(it)->r(im+) ----------
c     -- input 4 parameters -----
170 continue
      if (rhs(it) .lt. 0.0) stop 9999
      dq=sqrt(rhs(it))
      frf=1.0
      do 120 ii=1,4
       i=it+1-ii
       fo(i)=frf*dexp(-dq*r(i))
       go(i)=-dq*fo(i)
120 continue
      i=it-3
       FF1=DEXP(X(I))*Go(I)
       FF2=DEXP(X(I+1))*Go(I+1)
       FF3=DEXP(X(I+2))*Go(I+2)
       FF4=DEXP(X(I+3))*Go(I+3)
       GG1=rhs(i)*Fo(I)
       GG2=rhs(i+1)*Fo(I+1)
       GG3=rhs(i+2)*Fo(I+2)
       GG4=rhs(i+3)*Fo(I+3)
c     ********** inward loop 200 ******
      do 200 ii=3, it-im-1
       i=it-ii
C       ------- PREDICT -------
        FPR=Fo(I)-dx/24.0*(55.0*FF1-59.0*FF2+37.0*FF3-9.0*FF4)
        GPR=Go(I)-dx/24.0*(55.0*GG1-59.0*GG2+37.0*GG3-9.0*GG4)
C       ------- CORRECT -------
        FF0=DEXP(X(I-1))*GPR
        GG0=rhs(i-1)*FPR
        FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
        GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
        DO 1111 J=1, 2
         FF0=DEXP(X(I-1))*GCR
         GG0=rhs(I-1)*FCR
         FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
         GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
1111   CONTINUE
        Fo(I-1)=FCR
        Go(I-1)=GCR
C       ------------ PARAMETER EXCHANGE FOR THE NEXT LOOP -------
        FF4=FF3
        FF3=FF2
        FF2=FF1
        FF1=DEXP(X(I-1))*Go(I-1)
        GG4=GG3
        GG3=GG2
        GG2=GG1
        GG1=rhs(i-1)*Fo(I-1)
200 CONTINUE
c     ------------- integral of f r=0 to classical t. p.--------------
      sum=0.0d0+r(1)*f(1)*f(1)/2.0d0
      do 220 i=2, im-1, 2
       sum=sum+dx/3.0d0*(dexp(x(i-1))*f(i-1)**2+4.0d0*dexp(x(i))
     1      *f(i)**2+dexp(x(i+1))*f(i+1)**2)
220 continue
c     --------- integral of fo r=im to infinity ------
      sumout=0.0+fo(it)*fo(it)/2.0/dq
      do 225 i=im+1, it-1, 2
       sumout=sumout+dx/3.0d0*(dexp(x(i-1))*fo(i-1)**2
     &        +4.0d0*dexp(x(i))*fo(i)**2+dexp(x(i+1))*fo(i+1)**2)
225 continue
c     ------- energy correction estimate see hatree's book ---
      dele=(-g(im)/f(im)+go(im)/fo(im))/(sum/f(im)/f(im)+sumout
     &       /fo(im)/fo(im))
      if (abs(dele/e(inls)) .gt. 0.4) dele=0.8*abs(e(inls)/dele)*dele
      e(inls)=e(inls)-dele/2.0
      write (6,*) 'iscf/inls,iflag', iscf,inls,iflag,e(inls),
     $            dele/e(inls)
      if (abs(dele/e(inls)) .ge. 1.0d-5) go to 2010
c     --- if energy conv. is not satisfied go to 2010 --------
c      ---- matching the f and fo at classical t. p. ----
301     frfold=frf
      if (abs((f(im)-fo(im))/f(im)) .le. 1.0d-6) go to 300
      frf = frfold*f(im)/fo(im)
c     --------- recal fo from calssical t. p. to infinity
      do 121 ii=1,4
       i=it+1-ii
       fo(i)=frf*dexp(-dq*r(i))
       go(i)=-dq*fo(i)
121 continue
      i=it-3
       FF1=DEXP(X(I))*Go(I)
       FF2=DEXP(X(I+1))*Go(I+1)
       FF3=DEXP(X(I+2))*Go(I+2)
       FF4=DEXP(X(I+3))*Go(I+3)
       GG1=rhs(i)*Fo(I)
       GG2=rhs(i+1)*Fo(I+1)
       GG3=rhs(i+2)*Fo(I+2)
       GG4=rhs(i+3)*Fo(I+3)
c     ********** inward loop 201 ******
      do 201 ii=3, it-im-1
       i=it-ii
C       ------- PREDICT -------
        FPR=Fo(I)-dx/24.0*(55.0*FF1-59.0*FF2+37.0*FF3-9.0*FF4)
        GPR=Go(I)-dx/24.0*(55.0*GG1-59.0*GG2+37.0*GG3-9.0*GG4)
C       ------- CORRECT -------
        FF0=DEXP(X(I-1))*GPR
        GG0=rhs(i-1)*FPR
        FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
        GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
        DO 1112 J=1, 2
         FF0=DEXP(X(I-1))*GCR
         GG0=rhs(I-1)*FCR
         FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
         GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
1112   CONTINUE
        Fo(I-1)=FCR
        Go(I-1)=GCR
C       ------------ PARAMETER EXCHANGE FOR THE NEXT LOOP -------
        FF4=FF3
        FF3=FF2
        FF2=FF1
        FF1=DEXP(X(I-1))*Go(I-1)
        GG4=GG3
        GG3=GG2
        GG2=GG1
        GG1=rhs(i-1)*Fo(I-1)
201 CONTINUE
       go to 301
c     ------------- integral of f*f r=0 to infinity--------------
300 continue
      sum=0.0d0+r(1)*f(1)*f(1)/2.0d0
      do 221 i=2, im-1, 2
       sum=sum+dx/3.0d0*(dexp(x(i-1))*f(i-1)**2+4.0d0*dexp(x(i))
     1      *f(i)**2+dexp(x(i+1))*f(i+1)**2)
221 continue
      sumout=0.0+fo(it)*fo(it)/2.0/dq
      do 226 i=im+1, it-1, 2
       sumout=sumout+dx/3.0d0*(dexp(x(i-1))*fo(i-1)**2
     &        +4.0d0*dexp(x(i))*fo(i)**2+dexp(x(i+1))*fo(i+1)**2)
226 continue
      sum=sum+sumout
c      ------------------- normalization -----------
      do i=1, im
       fw(inls,i)=f(i)/sqrt(sum)
      enddo
      do i=im+1,it
       fw(inls,i)=fo(i)/sqrt(sum)
      enddo
      if (it .eq. 651) go to 227
       do i=it+1, 651
        fw(inls,i)=0.0
       enddo
227 continue
c     --------------- nomalization again -----------------
229 sum=0.0+r(1)*fw(inls,1)*fw(inls,1)/2.0d0
      do i=2, 650, 2
       sum=sum+dx/3.0d0*(dexp(x(i-1))*fw(inls,i-1)**2+
     & 4.0d0*dexp(x(i))*fw(inls,i)**2+dexp(x(i+1))*fw(inls,i+1)**2)
      enddo
      if (abs(sum-1.0d0) .le. 1.0d-10) go to 228
      do i=1, 651
       if (abs(fw(inls,i)) .gt. 1.0d-11)
     &      fw(inls,i)=fw(inls,i)/sqrt(sum)
      enddo
      go to 229
228 continue
c     ------------ summation of electron density ------
c     ------------- rosum = 4*pai*ro(r) ------------------
      do i=1, it
       rosum(nspin(inls),i)=rosum(nspin(inls),i)
     &         +wei(inls)*fw(inls,i)*fw(inls,i)/r(i)/r(i)
      enddo
2000 continue
c     ------ end of inls 2000 loop -------
c     ------ rosum check integ[r*r*rosum]=integ[4*pai*r*r*ro(r)]-----
       SEKI1=0.0+r(1)*r(1)*r(1)*rosum(1,1)/2.0
       seki2=0.0+r(1)*r(1)*r(1)*rosum(2,1)/2.0
       scrnc=seki1+seki2
       seki3=0.0+scrnc
       do i=1, 651
        ro(i)=rosum(1,i)+rosum(2,i)
       enddo
       DO 505 J=2, 650, 2
        SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*ROSUM(1,J-1)+4.0*
     &    DEXP(3.0*X(J))*ROSUM(1,J)+DEXP(3.0*X(J+1))*ROSUM(1,J+1))
        SEKI2=SEKI2+dx/3.0*(DEXP(3.0*X(J-1))*ROSUM(2,J-1)+4.0*
     &    DEXP(3.0*X(J))*ROSUM(2,J)+DEXP(3.0*X(J+1))*ROSUM(2,J+1))
        SEKI3=SEKI3+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*
     &    DEXP(3.0*X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
505 CONTINUE
       WRITE (6,*) 'atomic number=', z, 'rosum spin up=',SEKI1
       write (6,*) 'rosum spin down=',seki2, 'rosum=',seki3
C      ################# MODIFYING THE POTENTIAL ################
C      ******** HARTREE POTENTIAL ****************
       PAI4=4.0*3.141592654
       DO 500 I=1, 651
        SEKI1=0.0D0
        SEKI2=0.0D0
        IF (MOD(I,2)) 520,520,600
C       ------------ FOR EVEN I ----------
520   IF (I .EQ. 2) GO TO 530
        DO 540 J=2, I-2, 2
         SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*DEXP(3.0*
     1         X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
540   CONTINUE
530   SEKI1=SEKI1+dx/2.0*(DEXP(3.0*X(I-1))*RO(I-1)+
     1        DEXP(3.0*X(I))*RO(I))
        SEKI1=SEKI1+SCRNC
        IF (I .EQ. 650) GO TO 560
        DO 550 J=I+1, 649, 2
         SEKI2=SEKI2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*DEXP(2.0*
     1         X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
550   CONTINUE
560   SEKI2=SEKI2+dx/2.0*(DEXP(2.0*X(650))*RO(650)+
     1        DEXP(2.0*X(651))*RO(651))
        SEKI2=R(I)*SEKI2
        GO TO 680
C       ------------- FOR ODD I -----------
600   IF (I. EQ. 1) GO TO 620
        DO 610 J=2, I-1, 2
         SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*DEXP(3.0*
     1         X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
610   CONTINUE
620   SEKI1=SEKI1+SCRNC
        IF (I .EQ. 651) GO TO 640
        DO 650 J=I+1, 650, 2
         SEKI2=SEKI2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*DEXP(2.0*
     1         X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
650   CONTINUE
640   SEKI2=R(I)*SEKI2
C       ******** spin-dependent EXCHANGE-CORRELATION POTENTIAL ******
c         hedin lundqvist form J. Phys. C, 4 2064
c         von barth & hedin form: J.Phys.c 1972 5 p1629
c         vosko wilk nusair Can J Phys 1980 58 p1200
        pai=acos(-1.0)
680   IF (RO(I)) 690, 690, 700
690   PXup=0.0
        pxdown=0.0
        PCup=0.0
        pcdown=0.0
        GO TO 751
c       ------- rosum=(4*pai)*ro(r) -----------------
700   RS=(3.0/RO(I))**(1.0/3.0)

c       ------------ spin polarization spol--------------
        if (abs(rosum(1,i)-rosum(2,i)) .le. 1.0d-12) then
         spol=0.0
         fspol=0.0
         dfspol=0.0
         go to 704
        endif
        spol=(rosum(1,i)-rosum(2,i))/ro(i)
        fspol=((1.0+spol)**(4.0/3.0)+(1.0-spol)**(4.0/3.0)-2.0)
     &        /(2.0**(4.0/3.0)-2.0)
        dfspol=(4.0/3.0)*((1.0+spol)**(1.0/3.0)-(1.0-spol)**(1.0/3.0))
     &         /(2.0**(4.0/3.0)-2.0)
c       ------------ exchange ----------------------------------
704   exunpo=-(3.0/4.0/pai)*(9.0*pai/4.0)**(1.0/3.0)/rs
        expo=exunpo*(2.0)**(1.0/3.0)
        pxunpo=(4.0/3.0)*exunpo
        pxpo=(4.0/3.0)*expo
c        ----- unpo:unpolalized ,po:polarized,ex:exchange enrgy
c        ----- px:exchange potential
        ex=exunpo+(expo-exunpo)*fspol
        pxup=pxunpo+fspol*(pxpo-pxunpo)+(expo-exunpo)*
     &       (1.0-spol)*dfspol
        pxdown=pxunpo+fspol*(pxpo-pxunpo)+(expo-exunpo)*
     &       (-1.0-spol)*dfspol
c       ------------- correlation --------------------
c       --- Moruzzi Janak Williams form:Hedin-Lundqvist ----
        ahlunpo=21.0d0
        chlunpo=0.0225d0
        xhl=rs/ahlunpo
        alnunpo=log(1.0d0+1.0d0/xhl)
        pcunpo=-chlunpo*alnunpo
        ecunpo=-chlunpo*( (1.0d0+xhl**3)*log(1.0d0+1.0d0/xhl)+
     &        xhl/2.0d0-xhl**2-1.0d0/3.0d0)

        ahlpo=ahlunpo/2.0d0
        chlpo=chlunpo*2.0d0**(4.0/3.0)
        xhl=rs/ahlpo
        alnpo=log(1.0d0+1.0d0/xhl)
        pcpo=-chlpo*alnpo
        ecpo=-chlpo*( (1.0d0+xhl**3)*log(1.0d0+1.0d0/xhl)+
     &        xhl/2.0d0-xhl**2-1.0d0/3.0d0)

ca        if (rs-1.0) 701,702,702
ca 701   pcunpo=0.0311*log(rs)+(-0.048-0.0311/3.0)+
ca     &     (2.0/3.0)*0.0020*rs*log(rs)+(2.0*(-0.0116)-0.0020)/3.0*rs
ca        pcpo=0.01555*log(rs)+(-0.0269-0.01555/3.0)+
ca     &     (2.0/3.0)*0.0007*rs*log(rs)+(2.0*(-0.0048)-0.0007)/3.0*rs
ca        ecunpo=0.0311*log(rs)-0.048+0.0020*rs*log(rs)-0.0116*rs
ca        ecpo=0.01555*log(rs)-0.0269+0.0007*rs*log(rs)-0.0048*rs
ca        go to 750
ca 702   ecunpo=-0.1423/(1.0+1.0529*sqrt(rs)+0.3334*rs)
ca        ecpo=-0.0843/(1.0+1.3981*sqrt(rs)+0.2611*rs)
ca        pcunpo=ecunpo*(1.0+(7.0/6.0)*1.0529*sqrt(rs)+
ca     &     (4.0/3.0)*0.3334*rs)/(1.0+1.0529*sqrt(rs)+0.3334*rs)
ca        pcpo=ecpo*(1.0+(7.0/6.0)*1.3981*sqrt(rs)+
ca     &     (4.0/3.0)*0.2611*rs)/(1.0+1.3981*sqrt(rs)+0.2611*rs)

750      pcup=pcunpo+fspol*(pcpo-pcunpo)+(ecpo-ecunpo)
     &       *(1.0-spol)*dfspol
          pcdown=pcunpo+fspol*(pcpo-pcunpo)+(ecpo-ecunpo)
     &       *(-1.0-spol)*dfspol
c       --------- spin-dependent potential r*v(r)sigma -----
751   rvsig(1,i)=-Z+SEKI1+SEKI2+R(I)*(PCup+PXup)
        rvsig(2,i)=-Z+SEKI1+SEKI2+R(I)*(PCdown+PXdown)
        P1(1,I)=rv(1,I)
        p1(2,i)=rv(2,i)
C       ----------- CONSTRUCT THE NEW POTENTIAL ---------
        do jj=1, 2
         DENOM=P3(jj,I)+rvsig(jj,i)-P1(jj,I)-P2(jj,I)
         IF (DENOM .LE. 1.0D-10) GO TO 800
         IF (P1(jj,I)-rvsig(jj,i) .LE. 1.0D-10) GO TO 800
         AJUD=((P3(jj,I)*rvsig(jj,i)-P1(jj,I)*P2(jj,I))/DENOM
     &      -rvsig(jj,i))/(P1(jj,I)-rvsig(jj,i))
         IF (AJUD .GE. 0.5) GO TO 800
         rv(jj,I)=(P3(jj,I)*rvsig(jj,i)-P1(jj,I)*P2(jj,I))/DENOM
         IF (AJUD .LE. 0.0) rv(jj,I)=rvsig(jj,i)
         GO TO 820
800    rv(jj,I)=(P1(jj,I)+rvsig(jj,i))/2.0
C       -------------- POTENTIAL PARAMETER-----------------
820    IF (DABS(rv(jj,I)-P1(jj,I)) .GE. PJUD)
     &         PJUD=DABS(rv(jj,I)-P1(jj,I))
        enddo
500 CONTINUE
C      ##########################################################
       DO 420 I=1, 651
        do j=1, 2
         P3(j,I)=P1(j,I)
         P2(j,I)=rvsig(j,i)
         rvsig(j,i)=0.0
        enddo
420 CONTINUE
C      -------- POTENTIAL CONVERGENCE ---------------------
       IF (PJUD. LE. 5.0D-5) GO TO 3000
C      -------- ESTIMATE OF ENERGY CHANGE BY PURTURBATION -----
       DO 900 I=1, NO
        IF (wei(I) .EQ. 0.0) GO TO 900
        FF0=0.0D0
        do j=1, 651
         f(j)=fw(i,j)
        enddo
        DO 920 J=2, 650, 2
         FF1=(F(J-1)**2.0)*(rv(nspin(i),J-1)-P1(nspin(i),J-1))
         FF2=(F(J)**2.0)*(rv(nspin(i),J)-P1(nspin(i),J))
         FF3=(F(J+1)**2.0)*(rv(nspin(i),J+1)-P1(nspin(i),J+1))
         FF0=FF0+dx/3.0*(FF1+4.0*FF2+FF3)
920   CONTINUE
        E(I)=E(I)+FF0/2.0
900 CONTINUE
C      --------- CLEAR THE SUMMING PARAMETERS ---------------
       DO 410 I=1, 651
        ROSUM(1,I)=0.0
        rosum(2,i)=0.0
410 CONTINUE
       SCRNC=0.0D0
       PJUD=0.0D0
c       OPEN(1, FILE='RSCFEN.DAT')
c       WRITE (1,3010) (E(I), I=1, NO)
c       CLOSE (1)
c       OPEN(1, FILE='RSCFPO.DAT')
c       WRITE (1,3010) (P(I), I=1, 651)
c       CLOSE (1)
1000 CONTINUE
C     ************* SCF LOOP END **********************************
c     --------- kinetic energy ---------------------------
3000 continue
      do 3002 inls=1, 40
       ek(inls)=0.0
       if (wei(inls) .le. 1.0d-4) go to 3002
       do 3004 i=1, 651
        rhs(i)=(e(inls)-rv(nspin(inls),i)/r(i))*fw(inls,i)*fw(inls,i)
3004 continue
       ek(inls)=0.0+r(1)*rhs(1)/2.0
       do i=2, 650, 2
        ek(inls)=ek(inls)+dx/3.0d0*(dexp(x(i-1))*rhs(i-1)+
     &          4.0d0*dexp(x(i))*rhs(i)+dexp(x(i+1))*rhs(i+1))
       enddo
3002 continue
C     ************** OUTPUT THE RESULTS CALCULATED ****************
      DO 3030 I=1, 5
       WRITE (6,3020)
3030 CONTINUE
3020 FORMAT (1H ,'!!!!!!!!!!!!! CALCULATION DONE !!!!!!!!!!!!!!!')
      OPEN(1, FILE='RSCFEN.DAT')
      WRITE (1,3010) (E(I), I=1, 40)
      WRITE (6,3010) (E(I), I=1, 40)
3010 FORMAT (5(1PD15.7))
      CLOSE (1)
      OPEN(1, FILE='RSCFPO.DAT')
      WRITE (1,3010) (rv(1,I), I=1, 651)
      write (1,3010) (rv(2,i), i=1, 651)
      CLOSE (1)
      OPEN(1, FILE='ROSUM.DAT')
      write (1,3010) (ro(i), i =1, 651)
      WRITE (1,3010) (ROSUM(1,I), I=1, 651)
      write (1,3010) (rosum(2,i), i=1,651)
      CLOSE (1)
      OPEN(1, FILE='KINETIC.DAT')
      WRITE (1,3010) (ek(i), i=1, 40)
      WRITE (6,*   ) PJUD
      CLOSE (1)
      open (1, file='wavef.dat')
       do inls=1, no
        write (1,3010) (fw(inls,i), i=1, 651)
       enddo
       close (1)
      stop
      end