Table of Contents

Muffin-Tin Potential

Updated! This code may be copied for editing and/or compilation. Just remember to remove any HTML web coding lines from the source before trying to compile. OR ...
Download the gzipped version - mufpot.v03.4.f.gz

Created: April 7, 1999 ---- Last Updated: Sept 18, 2002
By Mark D. Pauli

__________________________________________________________________________



C     This program was used 7/99 for mufpot output for phase shifts
C     Don't trust muf9.out.  It is inconsistent between different platforms
C     and does not always give the same real ps as muf8.out.  M. Pauli 10/27/98
C--------------------------------------------------------------------------
C  MAIN PROGRAM FOR COMPUTATION OF MUFFIN-TIN POTENTIAL
C  AFTER MATTHEISS.   REF - LOUCKS, 'APW METHOD', 1967
C  INCLUDING% HARTREE POTENTIAL, SLATER-TYPE STATISTICAL
C  EXCHANGE TERM, AND MADELUNG CORRECTION FOR IONIC
C  MATERIALS
C  DIMENSIONED FOR UP TO 220 ATOMS IN THE UNIT CELL, OF 4
C  INEQUIVALENT TYPES
C
C       TO CHANGE THE DIMENSIONING:
C              FOR INEQUIVALENT TYPES, CHANGE 4'S BELOW
C                     (ORIGINAL VALUE=2)
C              FOR NUMBER OF ATOMS IN UNIT CELL, CHANGE 220'S
C                     (ORIGINAL VALUE=12)
C              FOR NUMBER OF SHELLS USED IN THE CALCULATION, CHANGE
C              THE 25'S BELOW, AND CHANGE THE DATA VALUE OF MC TO
C              AGREE
C                     (ORIGINAL VALUE=15)
C
C              DO THESE CHANGES IN MAIN
C                                   DIMENSION STATEMENTS
C                                   FORMAT STATEMENTS 911,912,913
C                                   COMMON BLOCKS
C                                HSIN, CLEMIN
C                                   COMMON BLOCKS
C                                RONUC
C                                   COMMON BLOCKS
C                                   DIMENSION STATEMENT
C            PHS(20),DELC(20) IN PHCALC SETS THE LIMIT LIM=LMAX+1 <= 20
C
C  PROGRAM HAS OPTIONS TO CALCULATE SCATTERING PHASE SHIFTS
C  USING SUBROUTINE PHCALC WITH OR WITHOUT AN L.S TERM
C  OR TO CALCULATE ATOMIC EXAFS MATRIX ELEMENTS USING
C  SUBROUTINE RONUC.

      PROGRAM MUFPOT
      DIMENSION SIG(260,4),RHO(260,4),VH(260,4),VS(260,4),
     &          VMAD(260,4),RX(260)
      DIMENSION RC(3,3),RK(3,220),ZM(220),Z(4),ZC(4),RMT(4),
     &          NRR(4),NCON(4),NX(4)
      DIMENSION RKK(3,220)
      DIMENSION IA(25,4),NA(25,4),AD(25,4)
      DIMENSION ATL(260),XM(260)
      DIMENSION RRHO(260),FNNS(260,4)
      CHARACTER*8 TITLE(9),NAME(4,4),WFN,WFN1,WFN2,WFN3,WFN4
      COMMON /WK/ WK1(441),WK2(441),PS(260)
      COMMON /WF/ WFC(260,25,4),WC(25),LC(25)
      COMMON /REL/ IREL,ALDS,VBCON
      DATA NGRID,MC/260,25/,WFN1,WFN2,WFN3,WFN4/
     &     'HERMAN-S','CLEMENTI','HEXINPUT','POTENTIA'/
      INDEX(X)=20.0*(ALOG(X)+8.8)+2.0

      PI = 4.0 * atan(1.0)
      HAR = 27.2113834
C
C     muf8.out gives either the EXAFS matrix elements or the stable
C       version of the real scattering phase shifts
C     muf9.out gives the fancy version of the real and imaginary
C       scattering phase shifts and seems to have a bug.
      OPEN(4,FILE='wf.dat',STATUS='OLD')
      OPEN(5,FILE='xt.dat',STATUS='OLD')
      OPEN(6,FILE='muf6.all.header.out',STATUS='UNKNOWN')
      OPEN(7,FILE='muf7.out',STATUS='UNKNOWN')
      OPEN(8,FILE='muf8.out',STATUS='NEW')
      OPEN(9,FILE='muf9.out',STATUS='UNKNOWN')
C
C     INITIALISATION OF LOUCKS' EXPONENTIAL MESH
      X = -8.8
      DO IX = 1,NGRID
        RX(IX)=EXP(X)
        X=X+0.05
      end do
C
C     INPUT OF CRYSTALLOGRAPHIC DATA
C      SPA = LATTICE CONSTANT IN A.U.
C      RC(I,J) = I'TH COORDINATE OF THE J'TH AXIS OF UNIT CELL, IN UNITS OF SPA
C      RK(I,J) = I'TH COORDINATE OF THE J'TH ATOM IN UNIT CELL, IN UNITS OF SPA
C      NR = NUMBER OF INEQUIVALENT ATOMS IN UNIT CELL
C     FOR AN ATOM OF TYPE IR
C      NRR(IR) = NUMBER IN UNIT CELL
C      Z(IR) = ATOMIC NUMBER
C      ZC(IR) = VALENCE CHARGE
C      RMT(IR) = MUFFIN-TIN RADIUS
C
C     IPX=0 IF PHASE SHIFTS ONLY REQUIRED
C        =1 IF X-RAY PROPERTIES REQUIRED
C     IREL=1 IF L.S TERM REQUIRED IN PHASE SHIFT CALCULATION
      READ(5,800)TITLE
      WRITE(6,900)TITLE
      READ(5,801)SPA,ANG
      READ(5,801)((RC(I,J),I=1,3),J=1,3)
      DO I = 1,3
        DO J = 1,3
          RC(I,J)=SPA*RC(I,J)
        end do
      end do
      READ(5,802)IPX,IREL
      IF((IPX.EQ.1) .AND. (IREL.EQ.1)) then
C       CALL TO RONUC (IPX=1), INCONSISTENT WITH USE OF L.S TERM (IREL=1)
        IREL=0
      end if
      READ(5,803)NR
      DO IR = 1,NR
        DO I = 1,NGRID
          VH(I,IR)=0.0
          VS(I,IR)=0.0
          VMAD(I,IR)=0.0
          SIG(I,IR)=0.0
          RHO(I,IR)=0.0
        end do 
      end do
      JJ=0
      ZZ=0.0
      DO IR = 1,NR
        READ(5,800)(NAME(I,IR),I=1,4)
        READ(5,804)NRR(IR),Z(IR),ZC(IR),RMT(IR)
        ZZ=ZZ+ABS(ZC(IR))
        N=NRR(IR)
        DO J = 1,N
          JJ=JJ+1
          ZM(JJ)=-ZC(IR)
          READ(5,801)RK(1,JJ),RK(2,JJ),RK(3,JJ)
          IF(ANG.NE.0.0) THEN
C           The coordinates of each atom were given in the form stated in the
C           books by R. W. G. Wyckoff, "Crystal Structures", Interscience,
C           (1965).  There, they are given as fractions of the unit vectors.
C           Rexpress the atom coordinates with respect to the Cartesian 
C           coordinate system to which the directions of the unit cell vectors,
C           RC, are referred.
            DO I = 1,3
              RKK(I,JJ)=RK(1,JJ)*RC(I,1) + RK(2,JJ)*RC(I,2) +
     &                  RK(3,JJ)*RC(I,3)
            end do
            DO I = 1,3
              RK(I,JJ)=RKK(I,JJ)
            end do
          ELSE
            DO I = 1,3
              RK(I,JJ)=RK(I,JJ)*SPA
            end do
          ENDIF
        end do
      end do
C 
C     N = TOTAL NUMBER OF ATOMS IN UNIT CELL
C     AV = TOTAL VOLUME OF UNIT CELL
C     OMA = ATOMIC VOLUME
C     RWS = WIGNER-SEITZ RADIUS
      N=JJ
      RCC1=RC(2,2)*RC(3,3)-RC(3,2)*RC(2,3)
      RCC2=RC(3,2)*RC(1,3)-RC(1,2)*RC(3,3)
      RCC3=RC(1,2)*RC(2,3)-RC(2,2)*RC(1,3)
      AV=ABS(RC(1,1)*RCC1+RC(2,1)*RCC2+RC(3,1)*RCC3)
      OMA=AV/FLOAT(N)
      RWS=(0.75*OMA/PI)**(1.0/3.0)
      WRITE(6,901)((RC(I,J),I=1,3),J=1,3)
      WRITE(6,902)AV,OMA,RWS
      JJ=0
      DO IR = 1,NR
        WRITE(6,903)IR,(NAME(I,IR),I=1,4),NRR(IR)
        INR=NRR(IR)
        DO IIR = 1,INR
          JJ=JJ+1
          WRITE(6,904)(RK(I,JJ),I=1,3)
        end do
        WRITE(6,905)Z(IR),ZC(IR),RMT(IR)
      end do
C     WRITE(6,906)(RX(IX),IX=1,NGRID)
C
C     FOR EACH ATOMIC TYPE, READ IN ATOMIC WAVEFUNCTIONS FOR NEUTRAL/positive
C     ATOM, IN HERMAN-SKILLMAN FORM (OUTPUT FROM PROGRAM 'HERSKL') OR
C     IN CLEMENTI FORM, OR CHARGE DENSITY FROM RELATIVISTIC SCF PROGRAM
C     'HEX', OR INPUT MUFFIN-TIN POTENTIAL DIRECTLY IN FINAL FORM.
C     RHO = 4*PI*CHARGE DENSITY * RADIUS**2
C     READ IN NHM THE MAXIMUM MUFFIN-TIN ZERO SAMPLING PARAMETER,
C     ALPHA FOR THE STATISTICAL EXCHANGE TERM(see Karlheinz Schwarz's
C     paper PRB V5,N7, 1/4/72, on Optimization of ALPHA),
C     ATDIAM OF ATOM (used in Im(E) for PHCALC),
C     FLOMT FOR THE FERMI LEVEL OFFSET FROM THE MUFFIN TIN ZERO (used 
C     in calc. of Im(E) for PHCALC),
C     EXPRES FOR THE EXPERIMENTAL RESOLUTION (used in Im(E)),
C     RMTZ FOR THE VALUE OF THE MUFFIN-TIN ZERO IF YOU DON'T DESIRE
C     TO HAVE IT CALCULATED BY THE PROGRAM. (IF NHM IS LESS THAN ZERO,
C     THIS NUMBER WILL BE USED)
      READ(5,804)NHM,ALPHA,ATDIAM,FLOMT,EXPRES,RMTZ
      EXPRES=EXPRES/HAR
      VINT=0.0
      MIX=0
      DO IR = 1,NR
        READ(4,800)WFN
C        write(*,*)wfn
C       OPTION 1)  HERMAN-SKILLMAN INPUT
        IF(WFN.EQ.WFN1)CALL HSIN(RHO(1,IR),RX,NX(IR),NGRID,IR)
C       OPTION 2)  CLEMENTI INPUT
        IF(WFN.EQ.WFN2)CALL CLEMIN(RHO(1,IR),RX,NX(IR),NGRID,IR)
C       OPTION 3)  HEX INPUT
        IF(WFN.EQ.WFN3)CALL HEXIN(RHO(1,IR),RX,NX(IR),NGRID)
        IF(WFN.EQ.WFN4) THEN 
C         OPTION 4) Muffin-Tin Potential in Hartrees
C         Read in R*V(R) into WK2 and radial grid into WK1 for the given 
C         number of points, NTAB.  The Muffin-tin Zero value relative to
C         vacuum is supplied as VINT in RYDBERGS.
C         NOTE:  With this option, only one atomic type is permitted.  The
C                loop over atom types is broken and control jumps to calculation
C                of phase shifts.
C         The final form of the total potential must be in AU*Rydbergs and of
C         the radii must be in AU.
C         Enough points of R*V(R) must be included to make it past the
C         muffin-tin radius.
C         CURRENTLY SET TO READ R IN AU AND R*V(R) IN AU*HARTREES
          READ(4,801)VINT
          WRITE(6,907)VINT
          READ(4,802)NTAB
          DO IX = 1,NTAB
            READ(4,805) WK1(IX),WK2(IX)
          end do
C         The following loop should be uncommented if WK1 and WK2 are given in
C         the units used by Herman-Skillman:
C           X = R / AMU and  U(r) = -r*V(r)_Ryd / (2*Z) = -r*V(r)_Har / Z
C         WK1 is converted from X to R in AU.
C         WK2 is converted from U to R*V(R) in AU*HARTREES.
C         DO IX = 1,NTAB
C           WK1(IX)=WK1(IX)*0.885341377/EXP(ALOG(Z(IR))/3.0)
C           WK2(IX)=-Z(IR)*WK2(IX)
C         end do
          NX(IR)=NGRID
          CALL CHGRID(WK2,WK1,NTAB,ATL,RX,NX(IR))
          NJC=NX(IR)
          DO JC = 1,NJC
C           Storing just the total potential in Rydbergs and adjusting the value
C           to be relative to the given muffin-tin zero.
            SIG(JC,IR)=2.0*ATL(JC)/RX(JC) - VINT
          end do
          GOTO 14
        end if

C       ******************************************************************
C       CALCULATE THE POTENTIAL FOR OPTIONS 1,2,3
C       ******************************************************************
C       RHO IS NORMALISED USING TOTAL ELECTRONIC CHARGE ON THE ATOM
C       CALCULATED BY THE TRAPEZOIDAL RULE
        NIX=NX(IR)
        MIX=MAX0(NIX,MIX)
        ZE=0.0
        W1=0.025*RHO(1,IR)*RX(1)
        DO IX = 2,NIX
          W2=0.025*RHO(IX,IR)*RX(IX)
          ZE=ZE+W1+W2
          W1=W2
        end do
        ZNT=Z(IR)-ZC(IR)
        ANORM=ZNT/ZE
        DO IX = 1,NIX
          RHO(IX,IR)=RHO(IX,IR)*ANORM
        end do
C       SOLVE POISSON'S EQUATION
C       SIG = COULOMB POTENTIAL
C       RHO = 4*PI*CHARGE DENSITY
        CALL POISSON(RHO(1,IR),ZNT,NIX,SIG(1,IR))
        X=-8.8
        DO IX = 1,NIX
          CE=EXP(-0.5*X)
          SIG(IX,IR)=CE*(-2.0*ZNT*CE+SIG(IX,IR))
          RHO(IX,IR)=RHO(IX,IR)*EXP(-2.0*X)
          X=X+0.05
        end do
        WRITE(6,908)(NAME(I,IR),I=1,4),ZE,ZNT,RX(NIX),NIX
C       WRITE(6,909)(SIG(IX,IR),IX=1,NIX)
      end do
C
C     DETAILS OF NEIGHBOURING SHELLS FOR EACH ATOMIC TYPE IR
C     NCON(IR) = NUMBER OF SHELLS INCLUDED
C     IA(J,IR) = ATOMIC TYPE IN J'TH SHELL
C     NA(J,IR) = NUMBER OF ATOMS IN J'TH SHELL
C     AD(J,IR) = DISTANCE TO J'TH SHELL
      RMAX=RX(MIX)
      CALL NBR(IA,NA,AD,NCON,NRR,NR,RC,RK,N,RMAX,MC)
      DO IR = 1,NR
        NC=NCON(IR)
        WRITE(6,910)(NAME(I,IR),I=1,4)
        IC=(NC-1)/8+1
        KC=0
        DO I = 1,IC
          JC=KC+1
          KC=MIN0(NC,KC+8)
          WRITE(6,911)(AD(J,IR),J=JC,KC)
          WRITE(6,912)(NA(J,IR),J=JC,KC)
          WRITE(6,913)(IA(J,IR),J=JC,KC)
        end do
      end do
C
C     CALCULATION OF THE MUFFIN-TIN POTENTIAL FOR EACH NEUTRAL
C     ATOM, FOLLOWING THE MATTHEISS PRESCRIPTION
      WRITE(6,914)NHM,ALPHA,RMTZ
      PD=6.0/(PI*PI)
      DO IR = 1,NR
C       SUMMING THE POTENTIALS FROM NEUTRAL ATOMS
C       VH = HARTREE POTENTIAL
        JRWS=INDEX(RMT(IR))
        CALL SUMAX(SIG,VH(1,IR),NX,NCON(IR),IA(1,IR),NA(1,IR),
     &             AD(1,IR),JRWS,NGRID,NR)
C       SUMMING THE CHARGE DENSITY ABOUT EACH ATOMIC TYPE
C       VS = TOTAL CHARGE DENSITY, THEN SLATER EXCHANGE TERM
        CALL SUMAX(RHO,VS(1,IR),NX,NCON(IR),IA(1,IR),NA(1,IR),
     &             AD(1,IR),JRWS,NGRID,NR)
        NIX=NX(IR)
        DO IX = 1,NIX
          VS(IX,IR)=-1.5*ALPHA*(PD*VS(IX,IR))**(1.0/3.0)
        end do
      end do
C
C     CALCULATE THE MUFFIN-TIN ZERO
      IF (NHM.GT.0) then
        CALL MTZ(NHM,SIG,RHO,RX,NGRID,RMT,NRR,NX,NR,RC,
     &           RK,N,VINT,ALPHA,AV,PI)
      else
C       AN ADDITION MADE BY G. HARP ON 24-DEC-89
C       SINCE CALCULATING THE MTZ IS THE LONGEST PART OF THE CALCULATION, ALLOW
C       FOR THE POSSIBILITY OF INPUTTING THIS PARAMETER FROM UNIT5.  IF NHM IS
C       NEGATIVE, THE PROGRAM USES THE VALUE READ FOR RMTZ, ABOVE. 
        VINT=RMTZ
      end if
C
C     THE MADELUNG CORRECTION FOR IONIC MATERIALS.   SUBROUTINE MAD
C     COMPUTES THE SPHERICALLY AND SPATIALLY AVERAGED FIELDS FOR
C     THE LATTICE OF POINT CHARGES ABOUT EACH ATOMIC TYPE
      IF(ZZ.NE.0)CALL MAD(VMAD,RX,NGRID,RMT,NRR,NX,NR,
     &                    RC,RK,ZM,N,AV,PI)
C
C     THE TOTAL MUFFIN-TIN POTENTIAL IS ACCUMULATED INTO SIG,
C     REFERRED TO THE MUFFIN-TIN ZERO IN RYDBERGS
      DO IR = 1,NR
        WRITE(6,915)(NAME(I,IR),I=1,4),RMT(IR)
        JRMT=INDEX(RMT(IR))
        DO IX = 1,JRMT-1
          SIG(IX,IR)=VH(IX,IR)+VS(IX,IR)+VMAD(IX,IR)-VINT
          WRITE(6,916)RX(IX),VH(IX,IR),VS(IX,IR),VMAD(IX,IR),SIG(IX,IR)
        end do
      end do

C     *******************************************************************
C     FOR ALL INPUT OPTIONS; WE HAVE THE TOTAL MUFFIN-TIN POTENTIAL, NOW
C     *******************************************************************
  14  DO 18 IR=1,NR
      READ(5,802)L,NE,E,DE
      IF ((L.LT.0) .AND. (IPX.EQ.0)) then
C       Skip phases for that atom type
      else
C       Calculate phases or matrix elements for atom type
C       Redefine the muffin-tin potential,SIG, into the form R*V(R) in
C       ATL in AU*Hartrees; and the muffin-tin zero, VINT, into WM in Hartrees.
        JRMT=INDEX(RMT(IR))
        ATL(1)=-Z(IR)
        XM(1)=0.0
        DO J = 2,JRMT
          XM(J)=RX(J-1)
          ATL(J)=SIG(J-1,IR)*XM(J)/2.0
        end do
        DO J = JRMT+1,NGRID
          XM(J)=RX(J-1)
          ATL(J)=0.0
        end do
        WM=VINT/2.0
        if (IPX .NE. 0) then
C         IF IPX=1 THE EXAFS ROUTINE RONUC IS TO BE CALLED
C         CALCULATE X-RAY SCATTERING FACTOR F(S)
C         S IS THE SCATTERING VECTOR FOR THE X-RAY
C           S IS 4.0*PI*SIN(THETA)/LAMBDA
C           F(S) IS 4.0*PI/S * (SINE FOURIER TRANSFORM OF R*CHARGE DENSITY)
C         NOTE THE DEFINITION OF RHO GIVEN ABOVE.
          SMIN=0.01
          SMAX=2.00
          SSTEP=SMIN
          IS1=1
          IS2=IFIX(SMAX/SSTEP+0.1)
          DO IX = 1,NIX
            RRHO(IX)=RHO(IX,IR)*RX(IX)/(4.0*PI)
          end do
          WRITE(6,917)(NAME(I,IR),I=1,4)
          DO I = 1,IS2
            S=SSTEP*I
            NS=3
            CALL FTSIN(RRHO,RX,NIX,S,FNNS(I,IR),NS)
            FNNS(I,IR)=FNNS(I,IR)*4.0*PI/S
            WRITE(6,918)S,FNNS(I,IR)
          end do
          IF (WFN.EQ.WFN1 .OR. WFN.EQ.WFN2) then
C           For HERMAN-Skillman or CLEMENTI inputs, calculate
C           Matrix Elements for X-ray absorption
            CALL RONUC(Z(IR),NGRID,ATL,XM,WM,RMT(IR),E,DE,NE,L,IR,PI)
          else 
C           Can not calculate matrix elements using HEXINPUT or POTENTIA input
            WRITE(6,919)
          end if
        else
C         IPX = 0, calculate scattering phase shifts
          WRITE(6,920)(NAME(I,IR),I=1,4)
          CALL PHCALC(Z(IR),NGRID,ATL,XM,WM,RMT(IR),E,DE,NE,L,
     &                ATDIAM,FLOMT,EXPRES,PI)
        end if
      end if
  18  CONTINUE
C
      STOP
C
 800  FORMAT(9A8)
 801  FORMAT(3F8.4)
 802  FORMAT(2I4,2F10.6)
 803  FORMAT(I4)
 804  FORMAT(I4,5F8.4)
 805  FORMAT(2(1pE20.12))
 900  FORMAT('MUFFIN-TIN POTENTIAL PROGRAM',//,9A8)
 901  FORMAT(//,' AXES OF UNIT CELL',/,(6X,3F8.4))
 902  FORMAT('   UNIT CELL VOLUME',F15.4,/,
     &       '   ATOMIC VOLUME',F18.4,/,
     &       '   WIGNER-SEITZ RADIUS',F12.4)
 903  FORMAT(//,' TYPE',I2,' ATOM  ::  ',4A8,/,
     &       I4,' ATOMS IN UNIT CELL')
 904  FORMAT(6X,3F8.4)
 905  FORMAT('   ATOMIC NUMBER',F18.4,/,
     &       '   VALENCE',F24.4,/,
     &       '   MUFFIN-TIN RADIUS',F14.4)
 906  FORMAT(//,' LOUCKS RADIAL MESH ',//,5(10F11.5,/))
 907  FORMAT('   MUFFIN-TIN ZERO',F16.5,' RYDBERGS')
 908  FORMAT(//,1X,4A8,/,' INTEGRATED ELECTRONIC CHARGE',F10.5,
     &       '  NORMALISED TO',F10.5,' ELECTRONS',/
     &       ' COULOMB POTENTIAL FOR ISOLATED NEUTRAL ATOM, '
     &       'OUT TO RADIUS',F10.5,', NX',I4,/)
 909  FORMAT(5(10E12.4,/))
 910  FORMAT(//,1X,4A8,/,' NEAREST NEIGHBOUR SHELLS ',/)
 911  FORMAT('   DISTANCE ',25(F8.4))
 912  FORMAT('   NUMBER   ',25(I5,3X))
 913  FORMAT('   TYPE     ',25(I5,3X))
 914  FORMAT(//,' MUFFIN-TIN ZERO CALCULATION',//,
     &          '       MAXIMUM MTZ MESH PARAMETER, NHM   ',I4,/,
     &          '   STATISTICAL EXCHANGE PARAMETER, ALPHA ',F8.4,/,
     &          '        MUFFIN-TIN ZERO PARAMETER, RMTZ  ',F8.4,
     &       ', RYDBERGS, READ FROM UNIT5')
 915  FORMAT(//,1X,4A8,/,' POTENTIALS IN RYDBERGS ',5X,
     &       ' OUT TO MUFFIN-TIN RADIUS ',F8.4,/,
     &       3X,' RADIUS',5X,'HARTREE POT',6X,'EXCHANGE',
     &       4X,'MADELUNG CORRCTN',4X,'TOTAL')
 916  FORMAT(F10.5,4(1pE16.5))
 917  FORMAT(//,1X,4A8,/,' EXAFS CALCULATION',//,
     &       '   X-RAY SCATTERING FACTOR',/,
     &       7X,'S',9X,'FNN(S)')
 918  FORMAT(3x,1pE10.4,3x,1pE11.5)
 919  FORMAT(//,' ****** INPUT INCONSISTENT WITH USE OF EXAFS ROUTINES',
     &        /,'   ****** CALL TO RONUC REJECTED. ',/)
 920  FORMAT(//,1X,4A8,/,' PHASE SHIFT CALCULATION',/)
C
      END




      SUBROUTINE HSIN(RHO,RX,NX,NGRID,IR)
      DIMENSION RHO(NGRID),RX(NGRID)
      REAL   NAME(4)
      COMMON /WK/ RR(441),RS(441),PS(260)
      COMMON /WF/ WFC(260,25,4),WC(25),LC(25)
C  ROUTINE FOR INPUT OF ATOMIC WAVEFUNCTIONS FROM HERMAN-SKILLMAN
C  TABLES, AND CALCULATION OF CHARGE DENSITY ON THE RADIAL MESH RX
C    RHO = 4*PI*SUM OVER STATES OF (MODULUS(WAVE FN)**2) *
C          RADIUS**2
C    NC = NUMBER OF ATOMIC STATES
C FOR EACH ATOMIC STATE I%
C    LC(I) = ANGULAR MOMENTUM
C    FRAC = FRACTIONAL OCCUPATION
C    WC(I) = NUMBER OF ELECTRONS
C    RR(IX)=RADIUS AT GRID POINT IX ON HERMAN-SKILLMAN GRID
C    WFC(IX,I,IR) = WAVEFUNCTION AT GRID POINT IX ON RADIAL GRID
      READ(4,140)NAME,Z,NC
      DO 1 IG=1,260
      RHO(IG)=0.0
      DO 1 IC=1,NC
   1  WFC(IG,IC,IR)=0.0
C INITIALISATION OF HERMAN-SKILLMAN radial MESH in AU R = X*AMU
      DR=0.00125*0.88534138/EXP(ALOG(Z)/3.0)
      RR(1)=0.0
      DO 2 I=2,441
      IF(MOD(I,40).EQ.2)DR=DR+DR
   2  RR(I)=RR(I-1)+DR
      NS=0
      DO 3 IC=1,NC
      READ(4,141)LC(IC),N,FRAC
      WC(IC)=FLOAT(2*LC(IC)+1)*FRAC*2.
      READ(4,142)(RS(IX),IX=1,N)
C  INTERPOLATION TO GRID RX
      NX=NGRID
      CALL CHGRID(RS,RR,N,PS,RX,NX)
      NS=MAX0(NS,NX)
      DO 3 IX=1,NX
      WFC(IX,IC,IR)=PS(IX)
C  CALCULATION OF CHARGE DENSITY
   3  RHO(IX)=RHO(IX)+WC(IC)*WFC(IX,IC,IR)**2
      WRITE(6,143)NAME
C     DO 4 IC=1,NC
C  4  WRITE(6,144)LC(IC),(WFC(IX,IC,IR),IX=1,NS)
C     WRITE(6,145)RX(NS),NS,(RHO(IX),IX=1,NS)
      RETURN
 140  FORMAT(4A8,/,F9.4,/,I4)
 141  FORMAT(I4,/,I4,/,F9.4)
 142  FORMAT(1P5E14.7)
 143  FORMAT(//,1X,4A8,/,' ATOMIC WAVEFUNCTIONS (HERMAN-SKILLMAN) ',
     &       'ON RADIAL MESH ',/)
 144  FORMAT(' L= ',I3,/,5(10E12.4,/))
 145  FORMAT(' CHARGE DENSITY OUT TO RADIUS ',F12.5,10X,
     &       ' NX',I4,/,5(10E12.4,/))
      END




      SUBROUTINE CLEMIN(RHO,RX,NX,NGRID,IR)
      DIMENSION RHO(NGRID),RX(NGRID)
      REAL   NAME(4)
      DIMENSION EX(20),FAC(20),FNT(20),NT(20)
      COMMON /WF/ WFC(260,25,4),WC(25),LC(25)
C  ROUTINE FOR INPUT OF WAVEFUNCTIONS IN THE CLEMENTI PARAMETRISED
C  FORM, AND CALCULATION OF CHARGE DENSITY ON RADIAL MESH RX
C    RHO = 4*PI*SUM OVER STATES OF (MODULUS(WAVE FN)**2) *
C          RADIUS**2
C    NC = NUMBER OF ATOMIC STATES
C  FOR EACH ATOMIC STATE I%
C    LC(I) = ANGULAR MOMENTUM
C    FRAC = FRACTIONAL OCCUPATION
C    WC(I) = NUMBER OF ELECTRONS
C    WFC(IX,I,IR) = WAVEFUNCTION AT GRID POINT IX
      READ(4,150)NAME
      READ(4,151)IPRINT
      READ(4,151)NC
      DO 1 IC=1,NC
      DO 1 IG=1,NGRID
   1  WFC(IG,IC,IR)=0.0
C  INPUT OF CLEMENTI PARAMETERS
      IC=1
   2  READ(4,151)NS
      IF(NS.LE.0)GOTO 8
      READ(4,152)(NT(I),EX(I),I=1,NS)
      DO 4 J=1,NS
      A=1.0
      B=2.0
      K=NT(J)
      C=FLOAT(K)
      KD=K+K
      DO 3 I=2,KD
      A=A*B
   3  B=B+1.0
   4  FNT(J)=EXP(-0.5*ALOG(A)+(C+0.5)*ALOG(2.0*EX(J)))
   5  READ(4,151)LC(IC)
      IF(LC(IC).LT.0)GOTO 2
      READ(4,153)(FAC(J),J=1,NS)
      READ(4,153)FRAC
      WC(IC)=2.0*FLOAT(2*LC(IC)+1)*FRAC
      DO 7 IX=1,NGRID
      A=0.0
      DO 6 K=1,NS
      EXX=EX(K)*RX(IX)
      IF(EXX.GT.80.0)GOTO 6
      A=A+FAC(K)*FNT(K)*(RX(IX)**(NT(K)))*EXP(-EXX)
   6  CONTINUE
   7  WFC(IX,IC,IR)=A
      IC=IC+1
      GOTO 5
C  CALCULATION OF CHARGE DENSITY
   8  DO 10 IX=1,NGRID
      SUM=0.0
      DO 9 IC=1,NC
   9  SUM=SUM+WC(IC)*WFC(IX,IC,IR)**2
      RHO(IX)=SUM
      IF(SUM.LT.1.0E-9)GOTO 11
  10  CONTINUE
  11  NX=IX
      IF(IPRINT.EQ.0)RETURN
      WRITE(6,154)NAME
C     DO 12 IC=1,NC
C 12  WRITE(6,155)LC(IC),(WFC(IX,IC,IR),IX=1,NGRID)
C     WRITE(6,156)RX(NX),NX,(RHO(IX),IX=1,NX)
      RETURN
 150  FORMAT(4A8)
 151  FORMAT(I4)
 152  FORMAT(I4,F11.5)
 153  FORMAT(5F11.5)
 154  FORMAT(1H1,4A8,32H ATOMIC WAVEFUNCTIONS (CLEMENTI))
 155  FORMAT(3H0L%,I3//5(10F11.5/))
 156  FORMAT(29H0CHARGE DENSITY OUT TO RADIUS,F12.5,10X,
     &       3HNX%,I4//5(10E12.4/))
      END




      SUBROUTINE HEXIN(RHO,RX,NX,NGRID)
      DIMENSION RHO(NGRID),RX(NGRID),RR(421),RS(421)
      REAL   NAME(4)
C  ROUTINE FOR INPUT OF WAVEFUNCTIONS CALCULATED FROM THE
C  RELATIVISTIC SELF-CONSISTENT FIELD PROGRAM HEX (MODIFIED
C  FROM THE CPC VERSION)
C
C  N             NUMBER OF MESH POINTS
C  J             NUMBER OF ATOMIC STATES
C  RN             MAXIMUM RADIUS OF MESH
C  H             INTERVAL OF THE LOGARITHMIC MESH. READ IN AS THE
C             RECIPROCAL
C  RR             INPUT RADIAL GRID
C  RS             RADIAL CHARGE DENSITY ON GRID RR
C  RX             LOUCK'S RADIAL GRID
C  RHO             CHARGE DENSITY INTERPOLATED ON RX
      READ(4,160)NAME,N,J,RN,H
      J=J+2
      DO 5 I=1,J
   5  READ(4,161)IDUMMY
      READ(4,162)(RS(I),I=1,N)
C  INITIALISATION OF INPUT MESH
      H=1.0/H
      D=EXP(H)
      R=RN/D**N
      DO 2 I=1,N
      R=R*D
   2  RR(I)=R
      NX=NGRID
      CALL CHGRID(RS,RR,N,RHO,RX,NX)
      WRITE(6,163)NAME,RX(NX),NX,(RHO(I),I=1,NX)
      RETURN
 160  FORMAT(4A8,/,2I5,F6.3,F4.0)
 161  FORMAT(1A4)
 162  FORMAT(1P5E14.7)
 163  FORMAT(1H ,4A8,6X,28HCHARGE DENSITY OUT TO RADIUS,F12.4,10X,
     &       3HNX%,I4,//,5(10E12.4/))
      END




      SUBROUTINE CHGRID(FX,X,NX,FY,Y,NY)
      DIMENSION FX(NX),X(NX),FY(NY),Y(NY)
C  PIECEWISE QUADRATIC INTERPOLATION FROM GRID X TO GRID Y,  BY
C  AITKEN'S DIVIDED DIFFERENCE SCHEME.   NX,NY ARE ARRAY DIMENSIONS%
C  NOTE THAT NY IS RESET IN CHGRID
      IY=1
      DO 2 IX=3,NX
   1  YY=Y(IY)
      IF(YY.GT.X(IX))GOTO 2
      A1=X(IX-2)-YY
      A2=X(IX-1)-YY
      A3=X(IX)-YY
      A12=(FX(IX-2)*A2-FX(IX-1)*A1)/(X(IX-1)-X(IX-2))
      A13=(FX(IX-2)*A3-FX(IX)*A1)/(X(IX)-X(IX-2))
      FY(IY)=(A12*A3-A13*A2)/(X(IX)-X(IX-1))
C CORRECTION-----------------------------------------------
      IY=IY+1
      IF(IY.GT.NY)GOTO 3
      GOTO 1
   2  CONTINUE
   3  NY=IY-1
      RETURN
      END




      SUBROUTINE POISSON(PSQ,Z,J,W)
      DIMENSION PSQ(J),W(J)
      DIMENSION EW(260),FW(260)
C CALCULATES COULOMB POTENTIAL FROM RADIAL DENSITY.
C TAKEN FROM LOUCKS' BOOK, APPENDIX 1
      A=1.0-0.0025/48.0
C EQ. A1.11
      B=-2.0-0.025/48.0
C EQ. A1.12
      EDL=EXP(0.025)
      C=0.0025/6.0
      C2=-B/A
      EW(1)=0.0
C EQ. A1.29
      FW(1)=EXP(0.025)
C EQ.A1.30
      X=-8.75
      ITOP=J-1
      DO 1 I=2,ITOP
      D=C*EXP(0.5*X)*(EDL*PSQ(I+1)+10.0*PSQ(I)+PSQ(I-1)/EDL)
C EQS. A1.13, A1.6
      FW(I)=C2-1.0/FW(I-1)
C EQ. A1.20
      EW(I)=(D/A+EW(I-1))/FW(I)
C EQ. A1.21
   1  X=X+0.05
      W(J)=2.0*Z*EXP(-0.5*X)
      DO 2 I=1,ITOP
      JC=J-I
   2  W(JC)=EW(JC)+W(JC+1)/FW(JC)
C EQ.A1.15
      RETURN
      END




      SUBROUTINE NBR(IA,NA,AD,NCON,NRR,NR,RC,RK,N,RMAX,MC)
      DIMENSION IA(MC,NR),NA(MC,NR),AD(MC,NR),NCON(NR),NRR(NR),
     &          RC(3,3),RK(3,N)
      DIMENSION RJ(3)
      RAD(A1,A2,A3)=SQRT(A1*A1+A2*A2+A3*A3)
C  ROUTINE TO SUPPLY NEAREST NEIGHBOUR DATA FOR ATOMS IN
C  A CRYSTAL STRUCTURE, GIVEN%
C  RC(I,J)% THE I'TH COORDINATE OF THE J'TH AXIS OF THE UNIT CELL
C  RK(I,J)% THE I'TH COORDINATE OF THE J'TH ATOM IN THE UNIT CELL
C  NRR(IR)% THE NUMBER OF TYPE-IR ATOMS IN THE UNIT CELL
C  THE INFORMATION RETURNED, FOR A TYPE-IR ATOM, IS
C  NCON(IR)% THE NUMBER OF NEAREST NEIGHBOUR SHELLS OF A TYPE-IR
C  ATOM INCLUDED, OUT TO A DISTANCE OF RMAX
C  IA(J,IR)% THE TYPE OF ATOMS IN THE J'TH NEIGHBOURING SHELL
C  NA(J,IR)% THE NUMBER OF ATOMS IN THE J'TH SHELL
C  AD(J,IR)% THE RADIUS OF THE J'TH SHELL
C  INITIALISATION
      RCMIN=1.0E6
      DO 1 I=1,3
   1  RCMIN=AMIN1(RCMIN,RAD(RC(1,I),RC(2,I),RC(3,I)))
      DO 2 IR=1,NR
      DO 2 IC=1,MC
      IA(IC,IR)=0
      NA(IC,IR)=0
   2  AD(IC,IR)=1.0E6
C  SEARCH OVER ADJACENT UNIT CELLS TO INCLUDE MC NEAREST NEIGHBOURS
      ITC=IFIX(RMAX/RCMIN)+1
      LIMC=ITC+ITC+1
      AS=-FLOAT(ITC+1)
      AX=AS
      DO 10 JX=1,LIMC
      AX=AX+1.0
      AY=AS
      DO 10 JY=1,LIMC
      AY=AY+1.0
      AZ=AS
      DO 10 JZ=1,LIMC
      AZ=AZ+1.0
      DO 3 J=1,3
   3  RJ(J)=AX*RC(J,1)+AY*RC(J,2)+AZ*RC(J,3)
C  RJ IS CURRENT UNIT CELL ORIGIN.  FOR EACH ATOM IN THIS UNIT CELL
C  FIND DISPLACEMENT R FROM KR-TYPE ATOM IN BASIC UNIT CELL
      J=0
      DO 10 JR=1,NR
      JNR=NRR(JR)
      DO 10 JJR=1,JNR
      J=J+1
      K=1
      DO 9 KR=1,NR
      R=RAD(RJ(1)+RK(1,J)-RK(1,K),RJ(2)+RK(2,J)-RK(2,K),
     &      RJ(3)+RK(3,J)-RK(3,K))
      IF(R.GT.RMAX)GOTO 9
C  COMPARE R WITH NEAREST NEIGHBOUR DISTANCES ALREADY FOUND
      IC=0
   4  IC=IC+1
      IF(IC.GT.MC)GOTO 9
      DR=R-AD(IC,KR)
C
C              HEY! THE FOLLOWING STATEMENT HAS BEEN REPLACED BY
C              G. HARP TO ALLOW A LITTLE MORE SLOP IN THE RADUIUS
C              CALCULATION: IF(ABS(DR).LT.1.0E-4)DR=0.0
C
      IF(ABS(DR).LT.1.0E-2)DR=0.0
      IF(DR)6,5,4
   5  IF(IA(IC,KR).NE.JR)GOTO 4
      NA(IC,KR)=NA(IC,KR)+1
      GOTO 9
   6  IF(IC.EQ.MC)GOTO 8
      IIC=IC+1
      DO 7 JJC=IIC,MC
      JC=MC+IIC-JJC
      IA(JC,KR)=IA(JC-1,KR)
      NA(JC,KR)=NA(JC-1,KR)
   7  AD(JC,KR)=AD(JC-1,KR)
   8  IA(IC,KR)=JR
      NA(IC,KR)=1
      AD(IC,KR)=R
   9  K=K+NRR(KR)
  10  CONTINUE
      DO IR = 1,NR
        NCON(IR)=0
        DO IC = 1,MC
          IF(NA(IC,IR).EQ.0)GOTO 12
          NCON(IR)=NCON(IR)+1
        end do
  12  end do
      RETURN
      END




      SUBROUTINE SUMAX(CHIAT,CHICRY,NX,NEIGHB,IA,NA,AD,JTOP,
     &                 NGRID,NR)
      DIMENSION CHIAT(NGRID,NR),CHICRY(NGRID),NX(NR),IA(1),NA(1),
     &          AD(1)
C  MODIFIED FROM LOUCKS APPENDIX 2% ROUTINE USES SPHERICAL EXPANSIONS
C  TO OBTAIN A TOTAL RADIAL DISTRIBUTION CHICRY FOR A PARTICULAR ATOM
C  FROM THE INDIVIDUAL RADIAL DISTRIBUTIONS CHIAT, FOR THAT ATOM AND
C  IT'S NEIGHBOURS.
C  NCON = NUMBER OF SHELLS INCLUDED ABOUT ATOM
C  IA(J) = ATOM TYPE IN J'TH SHELL ABOUT ATOM
C  NA(J) = NUMBER OF ATOMS IN J'TH SHELL ABOUT ATOM
C  AD(J) = DISTANCE OF J'TH SHELL FROM ATOM
      ICON=IA(1)
      JJCHI=NX(ICON)
      XJJCHI=JJCHI
      TOPX=0.05*(XJJCHI-1.)-8.8
      ICON=IA(1)
      DO 1 J=1,JTOP
   1  CHICRY(J)=CHIAT(J,ICON)
      JA=2
   2  X=-8.8
      I=1
      ICON=IA(JA)
   3  XINT=0.0
      ET=EXP(X)
      BLX=ALOG(AD(JA)-ET)
      IF (BLX-TOPX) 4,12,12
   4  JBL=2.0+20.0*(BLX+8.8)
      XJBL=JBL
      XBL=0.05*(XJBL-1.)-8.8
      G=XBL-BLX
      XINT=XINT+0.5*G*CHIAT(JBL,ICON)*(2.0-20.0*G)*EXP(2.*XBL)
      XINT=XINT+10.0*G*G*CHIAT(JBL-1,ICON)*EXP(2.0*(XBL-0.05))
      TLX=ALOG(AD(JA)+ET)
      IF (TLX-TOPX) 6,5,5
   5  JTL=NX(ICON)
      GO TO 9
   6  JTL=1.+20.0*(TLX+8.8)
      IF (JTL-JBL) 7,8,8
   7  FZN=CHIAT(JTL,ICON)*EXP(2.0*(XBL-0.05))
      FZ3=CHIAT(JBL,ICON)*EXP(2.0*XBL)
      FZ2=FZN+20.0*(FZ3-FZN)*(TLX-XBL+0.05)
      FZ1=FZN+20.0*(FZ3-FZN)*(BLX-XBL+0.05)
      XINT=0.5*(FZ1+FZ2)*(TLX-BLX)
      GO TO 12
   8  XJTL=JTL
      XTL=0.05*(XJTL-1.)-8.8
      C=TLX-XTL
      XINT=XINT+0.5*C*CHIAT(JTL,ICON)*(2.0-20.0*C)*EXP(2.0*XTL)
      XINT=XINT+10.0*C*C*CHIAT(JTL+1,ICON)*EXP(2.*(XTL+0.05))
   9  IF (JTL-JBL) 12,12,10
  10  XINT=XINT+0.025*CHIAT(JBL,ICON)*EXP(2.*XBL)
      XINT=XINT+0.025*CHIAT(JBL+1,ICON)*EXP(2.0*(XBL+0.05))
      JBL=JBL+1
      IF (JBL-JTL) 11,12,12
  11  XBL=XBL+0.05
      GO TO 10
  12  CHICRY(I)=CHICRY(I)+0.5*XINT*FLOAT(NA(JA))/
     &          (AD(JA)*ET)
      IF (I-JTOP) 13,14,14
  13  I=I+1
      X=X+0.05
      GO TO 3
  14  IF(JA-NEIGHB) 15,16,16
  15  JA=JA+1
      GO TO 2
  16  RETURN
      END




      SUBROUTINE MTZ(NHM,SIG,RHO,RX,NGRID,RMT,NRR,NX,NR,
     &               RC,RK,N,VINT,ALPHA,AV,PI)
      DIMENSION SIG(NGRID,NR),RHO(NGRID,NR),RX(NGRID),RMT(NR),
     &          NRR(NR),NX(NR),RC(3,3),RK(3,N)
      DIMENSION X(3),RB(3)
      RAD(A1,A2,A3)=SQRT(A1*A1+A2*A2+A3*A3)
      INDEX(XX)=20.0*(ALOG(XX)+8.8)+1.0
C  SUBROUTINE FOR CALCULATION OF THE MUFFIN-TIN ZERO LEVEL%
C  THE AVERAGE VALUE OF THE POTENTIAL BETWEEN THE MUFFIN-TIN
C  SPHERES IN THE UNIT CELL
C
      WRITE(6,170)
C  VOLT IS THE TRUE VOLUME OF THE REGION BETWEEN MUFFIN-TIN
C  SPHERES IN THE UNIT CELL
      VM=0.0
      DO 8 IR=1,NR
   8  VM=VM+FLOAT(NRR(IR))*RMT(IR)**3
      VOLT=AV-4.0*PI*VM/3.0
      PD=6.0/PI/PI
      VINT=0.0
      VHAR=0.0
      VEX=0.0
      NPOINT=0
      NINT=0
      DH=1.0
      NH=1
   1  AH=DH/2.0
      AX=-AH
      DO 7 IX=1,NH
      AX=AX+DH
      AY=-AH
      DO 7 IY=1,NH
      AY=AY+DH
      AZ=-AH
      DO 7 IZ=1,NH
      AZ=AZ+DH
      DO 2 I=1,3
   2  X(I)=AX*RC(I,1)+AY*RC(I,2)+AZ*RC(I,3)
      NPOINT=NPOINT+1
C  GIVES SAMPLE POINT X INSIDE THE UNIT CELL - TEST WHETHER
C  INTERSTITIAL
      DO 4 JX=1,2
      BX=FLOAT(JX-1)
      DO 4 JY=1,2
      BY=FLOAT(JY-1)
      DO 4 JZ=1,2
      BZ=FLOAT(JZ-1)
      DO 3 I=1,3
   3  RB(I)=BX*RC(I,1)+BY*RC(I,2)+BZ*RC(I,3)
      I=0
      DO 4 IR=1,NR
      INR=NRR(IR)
      DO 4 IIR=1,INR
      I=I+1
      XR=RAD(X(1)-RK(1,I)-RB(1),X(2)-RK(2,I)-RB(2),
     &       X(3)-RK(3,I)-RB(3))
      IF(XR.LT.RMT(IR))GOTO 7
   4  CONTINUE
C  WE HAVE AN INTERSTITIAL POINT
      NINT=NINT+1
C  SUM COULOMB AND EXCHANGE ENERGIES FROM ATOMS WITHIN 2 UNIT
C  CELLS AROUND THIS POINT
      SUMC=0.0
      SUME=0.0
      DO 6 JX=1,5
      BX=FLOAT(JX-3)
      DO 6 JY=1,5
      BY=FLOAT(JY-3)
      DO 6 JZ=1,5
      BZ=FLOAT(JZ-3)
      DO 5 I=1,3
   5  RB(I)=BX*RC(I,1)+BY*RC(I,2)+BZ*RC(I,3)
      J=0
      DO 6 JR=1,NR
      JNR=NRR(JR)
      DO 6 JJR=1,JNR
      J=J+1
      XR=RAD(RB(1)+RK(1,J)-X(1),RB(2)+RK(2,J)-X(2),RB(3)+RK(3,J)-X(3))
      J2=INDEX(XR)
      IF(J2.GE.NX(JR))GOTO 6
      J1=J2-1
      J3=J2+1
      X1=RX(J1)
      X2=RX(J2)
      X3=RX(J3)
      TERMC=(XR-X2)*(XR-X3)/(X1-X2)/(X1-X3)*SIG(J1,JR)
     &          +(XR-X1)*(XR-X3)/(X2-X1)/(X2-X3)*SIG(J2,JR)
     &          +(XR-X2)*(XR-X1)/(X3-X2)/(X3-X1)*SIG(J3,JR)
      TERME=(XR-X2)*(XR-X3)/(X1-X2)/(X1-X3)*RHO(J1,JR)
     &          +(XR-X1)*(XR-X3)/(X2-X1)/(X2-X3)*RHO(J2,JR)
     &          +(XR-X2)*(XR-X1)/(X3-X2)/(X3-X1)*RHO(J3,JR)
      SUMC=SUMC+TERMC
      SUME=SUME+TERME
   6  CONTINUE
C  CALCULATES HARTREE AND EXCHANGE CONTRIBUTIONS TO AVERAGE POTENTIAL
C  BETWEEN MUFFIN TINS
C
      VHAR=VHAR+SUMC
      VEX=VEX-1.5*ALPHA*(PD*SUME)**(1.0/3.0)
   7  CONTINUE
      DH=DH/2.0
      NH=NH+NH
      IF(NINT.EQ.0)GOTO 1
C  THE CURRENT MONTE-CARLO VOLUME FOR THE INTERSTITIAL REGION
C  IS VOLC
      VOLC=FLOAT(NINT)/FLOAT(NPOINT)*AV
      VHART=VHAR/FLOAT(NINT)
      VEXT=VEX/FLOAT(NINT)
      VINT=VHART+VEXT
      WRITE(6,171)VOLC,VINT
      IF(NH.LE.NHM)GOTO 1
C
      WRITE(6,172)NPOINT,VOLT,VOLC
      WRITE(6,173)VHART,VEXT,VINT
C
      RETURN
 170  FORMAT(/,'   MONTE CARLO ESTIMATES',/,
     &       7X,'VOLUME',5X,'MUF-TIN ZERO')
 171  FORMAT(F14.4,1X,F12.4)
 172  FORMAT('   SAMPLING WITH ',I9,' POINTS',//,
     &       '   INTERSTITAL REGION TRUE VOLUME ',F10.4,/,
     &       '               MONTE-CARLO VOLUME ',F10.4)
 173  FORMAT('    AVERAGE HARTREE POTENTIAL ',F9.5,' RYDBERGS',/,
     &       '   AVERAGE EXCHANGE POTENTIAL ',F9.5,' RYDBERGS',/,
     &       '    MUFFIN-TIN ZERO POTENTIAL ',F9.5,' RYDBERGS')
      END




      SUBROUTINE MAD(VMAD,RX,NGRID,RMT,NRR,NX,NR,RC,RK,ZM,N,AV,PI)
      DIMENSION VMAD(NGRID,NR),RX(NGRID),RC(3,3),RK(3,N),ZM(N),
     &          RMT(NR),NRR(NR),NX(NR)
      DIMENSION G(3,3),VMM(5),FR(5),RA(3),GA(3)
      DATA TEST/1.0E-4/
      RAD(A1,A2,A3)=SQRT(A1*A1+A2*A2+A3*A3)
C  SUBROUTINE MAD CALCULATES THE SPHERICALLY AND SPATIALLY AVERAGED
C  FIELDS FROM A LATTICE OF POINT CHARGES, AND TABULATES THEM ON
C  THE MESH RX
C  ** NB THIS ROUTINE WORKS IN HARTREES, BUT CONVERTS TO RYDBERGS **
C  RC(I,J) = THE I'TH COORDINATE OF THE J'TH AXIS OF THE UNIT CELL
C  RK(I,J) = THE I'TH COORDINATE OF THE J'TH ATOM IN THE UNIT CELL
C  VMAD(J,IR) = THE J'TH TABULATED VALUE OF THE SPHERICALLY AVERAGED
C  POTENTIAL ABOUT A TYPE-IR ATOM
C  ZM(K)=CHARGE ON THE K'TH ATOM
C  RMT(IR) = MUFFIN-TIN RADIUS OF A TYPE-IR ATOM
C  NR = NUMBER OF INEQUIVALENT ATOMS IN THE CELL
C  AV = VOLUME OF UNIT CELL
C  G(I,J) = I'TH COORDINATE OF THE J'TH RECIPROCAL LATTICE VECTOR
C  VMM(IR) = THE INTEGRAL OF THE POTENTIAL ABOUT A TYPE-IR ATOM
C  OUT TO THE MUFFIN-TIN RADIUS
C
      DO 1 IR=1,NR
      FR(IR)=0.0
      DO 1 J=1,NGRID
   1  VMAD(J,IR)=0.0
C  THE RECIPROCAL LATTICE IS DEFINED BY THREE VECTORS, G
      ATV=2.0*PI/AV
      G(1,1)=(RC(2,1)*RC(3,2)-RC(3,1)*RC(2,2))*ATV
      G(2,1)=(RC(3,1)*RC(1,2)-RC(1,1)*RC(3,2))*ATV
      G(3,1)=(RC(1,1)*RC(2,2)-RC(2,1)*RC(1,2))*ATV
C
      G(1,2)=(RC(2,2)*RC(3,3)-RC(3,2)*RC(2,3))*ATV
      G(2,2)=(RC(3,2)*RC(1,3)-RC(1,2)*RC(3,3))*ATV
      G(3,2)=(RC(1,2)*RC(2,3)-RC(2,2)*RC(1,3))*ATV
C
      G(1,3)=(RC(2,3)*RC(3,1)-RC(3,3)*RC(2,1))*ATV
      G(2,3)=(RC(3,3)*RC(1,1)-RC(1,3)*RC(3,1))*ATV
      G(3,3)=(RC(1,3)*RC(2,1)-RC(2,3)*RC(1,1))*ATV
C
C  MAXIMUM VALUE OF RK, AND MINIMUM VALUES OF RC,G - PRIOR TO
C  CHOOSING THE SEPARATION CONSTANT AL AND LIMITS FOR
C  SUMMATIONS
      RKMAX=0.0
      DO 2 J=1,N
   2  RKMAX=AMAX1(RKMAX,RAD(RK(1,J),RK(2,J),RK(3,J)))
      RCMIN=1.0E6
      GMIN=1.0E6
      DO 3 J=1,3
      RCMIN=AMIN1(RCMIN,RAD(RC(1,J),RC(2,J),RC(3,J)))
   3  GMIN=AMIN1(GMIN,RAD(G(1,J),G(2,J),G(3,J)))
C  AL IS CHOSEN TO GIVE EQUAL NUMBERS OF TERMS IN REAL AND
C  RECIPROCAL SPACE SUMMATIONS
      FAC1=TEST*ALOG(TEST)**4
      FAC2=(4.0*PI*RCMIN**4)/(AV*GMIN**4)
      AL=EXP(ALOG(FAC1/FAC2)/6.0)
      ITR=1+IFIX((AL*RKMAX-ALOG(TEST))/(AL*RCMIN))
      LIMR=ITR+ITR+1
      FAC1=4.0*PI*AL*AL/(AV*GMIN**4)
      ITG=1+IFIX(EXP(ALOG(FAC1/TEST)/4.0))
      LIMG=ITG+ITG+1
      WRITE(6,180)((G(I,J),I=1,3),J=1,3)
      WRITE(6,181)RCMIN,GMIN,RKMAX,TEST,AL
C
C  REAL SPACE SUMMATION
      WRITE(6,182)ITR
C  THE PREFACTORS FR FROM THE REAL SPACE SUMMATION ARE CALCULATED
      AS=-FLOAT(ITR)-1.0
      AX=AS
      DO 5 JX=1,LIMR
      AX=AX+1.0
      AY=AS
      DO 5 JY=1,LIMR
      AY=AY+1.0
      AZ=AS
      DO 5 JZ=1,LIMR
      AZ=AZ+1.0
      DO 4 I=1,3
   4  RA(I)=AX*RC(I,1)+AY*RC(I,2)+AZ*RC(I,3)
      DO 5 J=1,N
      K=1
      DO 5 KR=1,NR
      R=RAD(RA(1)+RK(1,J)-RK(1,K),RA(2)+RK(2,J)-RK(2,K),
     &      RA(3)+RK(3,J)-RK(3,K))
      IF(R.LT.1.0E-4)GOTO 5
      FR(KR)=FR(KR)+ZM(J)*EXP(-AL*R)/R
   5  K=K+NRR(KR)
      K=1
      DO 7 KR=1,NR
      X=RMT(KR)
      A=EXP(-AL*X)
      AI1=((1.0-A)/AL-X*A)/AL
      AI2=(X*0.5*(1.0/A+A)-0.5*(1.0/A-A)/AL)/AL/AL
      VMM(KR)=4.0*PI*(ZM(K)*AI1+AI2*FR(KR))
      NIX=NX(KR)
      DO 6 J=1,NIX
      X=RX(J)
      A=EXP(AL*X)
   6  VMAD(J,KR)=FR(KR)*0.5*(A-1.0/A)/(AL*X)+ZM(K)/(A*X)
   7  K=K+NRR(KR)
      WRITE(6,183)(VMM(KR),KR=1,NR)
C
C  NEXT COMES THE SUMMATION IN RECIPROCAL SPACE
      WRITE(6,184)ITG
      AS=-FLOAT(ITG)-1.0
      AX=AS
      DO 13 JX=1,LIMG
      AX=AX+1.0
      AY=AS
      DO 13 JY=1,LIMG
      AY=AY+1.0
      AZ=AS
      DO 13 JZ=1,LIMG
      AZ=AZ+1.0
      DO 8 I=1,3
   8  GA(I)=AX*G(I,1)+AY*G(I,2)+AZ*G(I,3)
      GM=RAD(GA(1),GA(2),GA(3))
      GS=GM*GM
      FAC1=0.0
      IF(GS.LT.1.0E-4)GOTO 13
      FAC1=4.0*PI*AL*AL/(AV*GS*(GS+AL*AL))
      K=1
      DO 12 KR=1,NR
      FAC2=0.0
      DO 10 J=1,N
      GR=0.0
      DO 9 I=1,3
   9  GR=GR+GA(I)*(RK(I,K)-RK(I,J))
  10  FAC2=FAC2+COS(GR)*ZM(J)
      X=RMT(KR)
      AI3=(SIN(GM*X)/GM-X*COS(GM*X))/GS
      VMM(KR)=VMM(KR)+4.0*PI*AI3*FAC1*FAC2
      NIX=NX(KR)
      DO I = 1,NIX
        X=RX(I)
        VMAD(I,KR)=VMAD(I,KR)+FAC1*FAC2*SIN(GM*X)/(GM*X)
      end do
  12  K=K+NRR(KR)
  13  CONTINUE
      WRITE(6,183)(VMM(KR),KR=1,IR)
C  REFER TO MUFFIN-TIN ZERO
      VM=0.0
      AMT=0.0
      DO 14 IR=1,NR
      VM=VM+FLOAT(NRR(IR))*RMT(IR)**3
  14  AMT=AMT+FLOAT(NRR(IR))*VMM(IR)
      AMT=AMT/(AV-4.0*PI*VM/3.0)
C  EXPRESS THE FINAL POTENTIAL IN RYDBERGS
      AMT=-2.0*AMT
      WRITE(6,185)AMT
      DO 15 KR=1,NR
      NIX=NX(KR)
      DO 15 J=1,NIX
  15  VMAD(J,KR)=2.0*VMAD(J,KR)-AMT
      RETURN
C
 180  FORMAT(//,' MADELUNG CORRECTION',//,
     &          '   RECIPROCAL LATTICE',/,(6X,3F8.4))
 181  FORMAT('   RCMIN ',F10.4,10X,'GMIN ',F10.4,10X,'RKMAX ',F10.4,
     &       10X,'TEST ',E12.4,/,' SEPARATION CONSTANT ',E12.4)
 182  FORMAT('   REAL SPACE SUMMATION ',11X,'ITR ',I3)
 183  FORMAT('   VMM (HARTREES) ',5E12.4)
 184  FORMAT('   RECIPROCAL SPACE SUMMATION ',5X,'ITG ',I3)
 185  FORMAT('   MADELUNG MUFFIN-TIN ZERO ',5E12.4)
C
      END




COMMENT THIS SUBROUTINE (SUPPORTED BY CP) FOURIER TRANSFORMS THE
C       INPUT FUNCTION FUN. FUN IS READ IN AT GRID POINTS GIVEN
C       BY THE ARRAY R,AND IS FOURIER TRANSFORMED WITH RESPECT TO
C       SIN(A*R). A IS A SIMPLE INPUT VARIABLE.
      SUBROUTINE FTSIN(FUN,R,NR,A,FTI,NS)
      DIMENSION F(5),P(5),D(5),FL(5),FR(5),CL(5,5),CR(5,5)
      DIMENSION FUN(260),R(260),SC(2)
      NSX=5
      NRX=260
      IF(NR.GT.NRX) GO TO 90
      NS=MIN0(NS,NSX)
      IR=0
      F(1)=FUN(1)
      P(1)=R(1)
      SC(1)=SIN(P(1)*A)/A
      SC(2)=COS(P(1)*A)/(A*A)
      FTI=F(1)*SC(2)
COMMENT. NS FIXES THE DEGREE OF THE POLYNOMIAL USED TO APPROXIMATE
C        FUN, AND HENCE THE NUMBER OF SUB-INTERVALS, WITHIN THE
C        INTERVAL OF THE R AXIS.
   30 IF(IR.GT.(NR-NS)) GO TO 70
      DO 31 JS=2,NS
      IJ=IR+JS
      P(JS)=R(IJ)
   31 F(JS)=FUN(IJ)
      CALL CP(P,D,CL,CR,NS,NSX)
      DO 41 JD=2,NS
      FL(JD)=0.0
      FR(JD)=0.0
      DO 41 JS=1,NS
      FL(JD)=FL(JD)+F(JS)*CL(JD,JS)
   41 FR(JD)=FR(JD)+F(JS)*CR(JD,JS)
C
   50 ASQ=A*A
      SLA=SC(1)
      SRA=SIN(P(NS)*A)/A
      SC(1)=SRA
      SLB=SC(2)
      SRB=COS(P(NS)*A)/ASQ
      SC(2)=SRB
      DO 51 JD=2,NS
      SSA=SLB
      SLB=-SLA/ASQ
      SLA=SSA
      SSA=SRB
      SRB=-SRA/ASQ
      SRA=SSA
   51 FTI=FTI-(FR(JD)*SRB-FL(JD)*SLB)
C
      F(1)=F(NS)
      P(1)=P(NS)
      IR=IR+NS-1
      GO TO 30
   70 FTI=(FTI-F(1)*SC(2))*A
      GOTO 91
   90 WRITE(6,*)'NUMBER OF R POINTS TOO LARGE'
   91 RETURN
      END




      SUBROUTINE CP(P,D,CL,CR,NS,M)
      DIMENSION P(NS),D(NS), CL(M,NS),CR(M,NS)
      GO TO (10,20,30,40), NS
   10 RETURN
   20 D(1)=1./(P(1)-P(2))
      D(2)=1./(P(2)-P(1))
      CL(2,1)=D(1)
      CL(2,2)=D(2)
      CR(2,1)=D(1)
      CR(2,2)=D(2)
      RETURN
   30 D(1)=1./((P(1)-P(2))*(P(1)-P(3)))
      D(2)=1./((P(2)-P(1))*(P(2)-P(3)))
      D(3)=1./((P(3)-P(1))*(P(3)-P(2)))
      CL(2,1)=(P(1)-P(2)+P(1)-P(3))*D(1)
      CL(2,2)=(P(1)-P(3))*D(2)
      CL(2,3)=(P(1)-P(2))*D(3)
      CR(2,1)=(P(3)-P(2))*D(1)
      CR(2,2)=(P(3)-P(1))*D(2)
      CR(2,3)=(P(3)-P(1)+P(3)-P(2))*D(3)
      CL(3,1)=2.*D(1)
      CL(3,2)=2.*D(2)
      CL(3,3)=2.*D(3)
      CR(3,1)=2.*D(1)
      CR(3,2)=2.*D(2)
      CR(3,3)=2.*D(3)
      RETURN
   40 D(1)=1./((P(1)-P(2))*(P(1)-P(3))*(P(1)-P(4)))
      D(2)=1./((P(2)-P(1))*(P(2)-P(3))*(P(2)-P(4)))
      D(3)=1./((P(3)-P(1))*(P(3)-P(2))*(P(3)-P(4)))
      D(4)=1./((P(4)-P(1))*(P(4)-P(2))*(P(4)-P(3)))
      CL(2,1)=((P(1)-P(3))*(P(1)-P(4))+(P(1)-P(2))*(P(1)-P(4))
     &                 +(P(1)-P(2))*(P(1)-P(3)))*D(1)
      CL(2,2)=(P(1)-P(3))*(P(1)-P(4))*D(2)
      CL(2,3)=(P(1)-P(2))*(P(1)-P(4))*D(3)
      CL(2,4)=(P(1)-P(2))*(P(1)-P(3))*D(4)
      CR(2,1)=(P(4)-P(2))*(P(4)-P(3))*D(1)
      CR(2,2)=(P(4)-P(1))*(P(4)-P(3))*D(2)
      CR(2,3)=(P(4)-P(1))*(P(4)-P(2))*D(3)
      CR(2,4)=((P(4)-P(2))*(P(4)-P(3))+(P(4)-P(1))*(P(4)-P(3))
     &                 +(P(4)-P(1))*(P(4)-P(2)))*D(4)
      CL(3,1)=2.*(P(1)-P(2)+P(1)-P(3)+P(1)-P(4))*D(1)
      CL(3,2)=2.*(P(1)-P(3)+P(1)-P(4))*D(2)
      CL(3,3)=2.*(P(1)-P(2)+P(1)-P(4))*D(3)
      CL(3,4)=2.*(P(1)-P(2)+P(1)-P(3))*D(4)
      CR(3,1)=2.*(P(4)-P(2)+P(4)-P(3))*D(1)
      CR(3,2)=2.*(P(4)-P(1)+P(4)-P(3))*D(2)
      CR(3,3)=2.*(P(4)-P(1)+P(4)-P(2))*D(3)
      CR(3,4)=2.*(P(4)-P(1)+P(4)-P(2)+P(4)-P(3))*D(4)
      CL(4,1)=6.*D(1)
      CL(4,2)=6.*D(2)
      CL(4,3)=6.*D(3)
      CL(4,4)=6.*D(4)
      CR(4,1)=6.*D(1)
      CR(4,2)=6.*D(2)
      CR(4,3)=6.*D(3)
      CR(4,4)=6.*D(4)
      RETURN
      END




      SUBROUTINE PHCALC(Z,NM,VM,XM,WM,RM,E1,DE,NE,LMAX,ATDIAM,FLOMT,
     &                  EXPRES,PI)
      COMPLEX DFI(200),FI(200),EE,DELC(20)
      DIMENSION STORE(200),PHS(20)
      DIMENSION XM(260),VM(260)
      COMMON /WS/KCF,JRCF,Y(200),CL,E,VL(200)
      COMMON /REL/IREL,ALDS,VBCON
      DATA ET/3.0/
COMMENT  ROUTINE TRANSFERS POTENTIAL TO H-F GRID
C        AND CALLS CALCFI TO CALCULATE PHASE SHIFT.
C  XM=LOUCK'S RADIAL GRID
C  VM=MUFFIN-TIN POTENTIAL X RADIAL DAISTANCE ON LOUCK'S GRID
C  Y=HERMAN-SKILLMAN GRID converted to radius (AU)
C  VL=POTENTIAL X RADIAL DISTANCE ON H-S GRID
C  L=ANGULAR MOMENTUM
C  E=ENERGY
C  PHS=PHASE SHIFT
C
C  PHASE SHIFTS OUTPUT ON UNIT 8
C  MUFFIN-TIN POTENTIAL ON H-S GRID OUPUT ON UNIT 7
      WRITE(6,190) Z,NM,WM,RM
      WRITE(6,191) (XM(J),VM(J),J=1,NM)
COMMENT. SET UP H-F GRID (HERMAN-SKILLMAN FORM).
      NMESH=200
      AMU=0.88534138/EXP(ALOG(Z)/3.0)
      DX=0.01*AMU
      ES=SQRT(2.0*ET)
      Y(1)=0.0
      J=2
      NF=NMESH
   2  KDL2=J+9
      DO 3 K=J,KDL2
   3  Y(K)=Y(K-1)+DX
      J=J+10
      IF(J-NMESH+10) 4,4,8
   4  IF(0.3/ES-DX) 6,6,5
   5  DX=DX+DX
      GO TO 2
   6  NF=J-1
      DO 7 K=J,NMESH
   7  Y(K)=Y(K-1)+DX
   8  CONTINUE
      CALL CHGRID(VM,XM,NM,VL,Y,NMESH)
      NMESH=200
      IXX=0
      DO 9 J=1,NMESH
      IF(Y(J).GT.RM .AND. IXX.EQ.0) THEN
      IXX=1
      KRMT=J
      ENDIF
   9  IF(Y(J).GT.RM) VL(J)=0.0000
      WRITE(7,192)NMESH,KRMT,NF
      WRITE(6,193)ET
      WRITE(6,194) (Y(J),VL(J),J=1,NMESH)
      WRITE(7,195) (Y(J),VL(J),J=1,NMESH)
      WRITE(6,196)
      LIM=LMAX+1
      S=0.5
      IF(IREL.NE.1)GOTO 30
      WRITE(6,197)
  30  E=E1
      DO 11 JE=1,NE
      DO 10 JL=1,LIM
      L=JL-1
      IF(IREL.NE.1)GOTO 40
      ALDS=FLOAT(L)*S
  40  CALL CALCFI(L,Z,STORE,DEL,NF,RM,PI)
789   format(F10.5,I5,F10.5)
C
C       SIGMA :: SUPPLIES THE SELF-ENERGY FOR ENERGY E; i.e. SIGMA = Im(E)
C       BY CONVENTION, NEGATIVE VALUES OF SIGMA CORRESPOND TO ABSORPTION.
C       Given u = A exp(ik.x)  with k = Re(k) + i Im(k)
C       Then I = |u|^2 = |A|^2 exp(-2 Im(k) x) where for Im(k) << Re(k)
C       Im(k) = Im(E - V) / sqrt(2 Re(E - V)).
C       By definition Im(k) = 1/lambda_amp = 1/(2 lambda_int) such that
C       exp(-2 Im(k) x) = exp(-2x/lambda_amp) = exp(-x/lambda_int).
C         ADDITION MADE BY G. HARP 23-JAN-90
C         MODIFIED BY M. PAULI 27-OCT-98
C
C       FIRST LAMBDA IS CALCULATED BY THE EMPIRICAL FORMULA OF SEAH AND DENCH 
C       (SOURCE OF THE "UNIVERSAL CURVE") Surf. and Interface Anal. 1,1 (1979).
C       i.e. LAMBDA  = 538/E^2 + 0.41*(aE)^0.5
C
C         WHERE     E = ELECTRON KINETIC ENERGY MEASURED FROM FERMI LEVEL (eV)
C                   a = ATOMIC DIAMETER (nm)
C              LAMBDA = ABSORPTION LENGTH FOR INTENSITY (monolayers)
C
C         SINCE THE ABOVE EQUATION MEASURES E FROM (A SOMEWHAT ARBITRARY)
C         "FERMI LEVEL", A FERMI LEVEL OFFSET FROM THE MUFFIN TIN ZERO (FLOMT)
C         IS INCLUDED AS A PARAMTER.  THIS MAINLY BECOMES IMPORTANT WHEN THE
C         ELECTRON KINETIC ENERGY IS SMALL.  (ONE MIGHT EXPECT THE FERMI LEVEL
C         TO TYPICALLY BE AROUND 5 eV BELOW VACUUM AND THE MUFFIN TIN CAN
C         BE ANYWHERE INSIDE THE VALENCE BAND REGION SOMEWHERE BELOW THE
C         FERMI LEVEL; SAY AROUND 8-10 eV or more BELOW VACUUM.  THEN A
C         TYPICAL VALUE FOR FLOMT IS AROUND 0.25 Ryd..)
C
C       SECOND, THE INITIAL SIGMA IS CALCULATED
C       i.e. SIGMA = Im(E) = SQRT(2*Re(E))/(2*LAMBDA)
C       THUS SIGMA = 0.5 / ( 538/HAR^2/SQRT(2)*E^(-2.5)*a
C                            + 0.41*SQRT(HAR*AU/10/2)*a^(1.5) )
C
C         WHERE   Re(E) = ELECTRON REAL KINETIC ENERGY (Har)
C                LAMBDA = ABSORPTION LENGTH FOR INTENSITY (bohr radius:AU)
C                   HAR = 27.211383 (eV)
C                    AU = 0.529177249 (Ang)
C
C       FINALLY, TO ACCOUNT FOR DETECTOR LOSS, THE EXPERIMENTAL RESOLUTION 
C       (EXPRES) IS QUADRATICALLY ADDED TO Im(E) TO GET AN EFFECTIVE SIGMA.
C       IF ATDIAM=0, WE WANT TO BYPASS ALL OF THIS AND JUST SET SIGMA=EXPRES;
C       (THIS FORCES SIGMA TO RETURN A CONSTANT STEP SIZE TO THE PROGRAM.)
C
C       THE PROGRAM DOES SOME UNIT CONVERSIONS.  INPUTS ARE SUCH THAT 
C         ATDIAM IS ENTERED IN Bohr Radius (AU)
C         FLOMT  IS ENTERED IN Hartrees
C         EXPRES IS ENTERED IN eV
C         THE GIVEN E IS THE ENERGY REFERENCE TO THE MUFFIN TIN ZERO.
C         CALCULATED SE IS ENERGY REFERENCE TO THE FERMI LEVEL.
C

        SE=E-FLOMT
       IF(ATDIAM.NE.0) THEN
          SIGMA=.5/(.51377*ATDIAM*(SE**(-2.5))+0.34789*(ATDIAM**(1.5)))
          SIGMA=SQRT(SIGMA**2+EXPRES**2)
       ELSE
          SIGMA=EXPRES
       ENDIF
       WRITE(6,198)SE,SIGMA
C
C
      EE=CMPLX(E,SIGMA)
      CALL CALPFI(FI,DFI,FI,DFI,DELC(JL),EE,L,1,0,0,Y,VL,NMESH,KRMT,Z,
     &            NF)
  10  PHS(JL)=DEL
C
C
C   ENERGIES AND PHASE-SHIFTS ARE WRITTEN FROM UNIT 9 ONTO
      WRITE(9,199) E,(REAL(DELC(JL)),JL=1,LIM)
      WRITE(9,199) E,(AIMAG(DELC(JL)),JL=1,LIM)
C   A FILE IN CORRECT FORMAT FOR INPUT TO THE XANES PROGRAM
C
C
      IF(S.LT.0.0)GOTO 13
      WRITE(6,200) E,(PHS(JL),JL=1,LIM)
      WRITE(8,199) E,(PHS(JL),JL=1,LIM)
      GOTO 11
  13  WRITE(6,201) E,(PHS(JL),JL=2,LIM)
      WRITE(8,202) E,(PHS(JL),JL=2,LIM)
  11  E=E+DE
      IF(IREL.NE.1.OR.S.LT.0.0)GOTO 50
      S=-0.5
      WRITE(6,203)
      GOTO 30
  50  RETURN
 190  FORMAT('   Z=',F9.4,'   NM=',I4,'   WM=',F8.4,'   RM=',F7.4,/)
 191  FORMAT('        XM          -XM.VM',/,200(/,1F12.8,1F14.5))
 192  FORMAT(3I5)
 193  FORMAT('  POTENTIAL (V(R)*R) IN HAR. ON H-S GRID FOR E=',F8.4/,
     &       '      ',3('R[AU]       VL          '),/)
 194  FORMAT(6E12.4)
 195  FORMAT(5E14.6)
 196  FORMAT(/,' PHASE SHIFTS',/)
 197  FORMAT(' SPIN UP'/)
 198  FORMAT('FOR E=',F8.4,'  , SIGMA=',F8.4)
 199  FORMAT(7(1pE12.5))
 200  FORMAT(' E=',1F12.5,10(/,5F12.5))
 201  FORMAT(' E=',1F12.5,/,4X,8(1H-),4F12.5,9(/,5F12.5))
 202  FORMAT(6(1pE12.5))
 203  FORMAT(/,' SPIN DOWN',/)
      END




      SUBROUTINE CALCFI(L,Z,STORE,DEL,NF,RHO,PI)
      DIMENSION STORE(200)
      DIMENSION A(2),B(2),C(2),D(2),F(2),RAD(3)
      COMMON /WS/K,JR,X(200),CL,E,ATP(200)
      COMMON /REL/IREL,ALDS,VBCON
      DATA CC/18779.688/
COMMENT  CALCULATES PHASE SHIFTS BY SEPARATING SCHROEDINGER
C        EQUATION INTO TWO FIRST ORDER EQUATIONS AND INTEGRATING
C        USING RUNGE-KUTTA-MERSON SUBROUTINE DA01A
C      L=ANGULAR MOMENTUM
C      E=ENERGY
C      DEL=PHASE SHIFT
     
      BIG=0.
      NMESH=200
      DO 29 J=1,NMESH
  29  STORE(J)=0.0
   1  CL=FLOAT(L)
      M=1
      SPLIT=1.0
   4  DX=X(2)-X(1)
COMMENT SPLIT IS THE DEGREE OF SUBDIVISION OF DX, THE BASIC
C       INTERVAL OF TABULATION. DEPENDING ON WHETHER THE OPTION
C       TO INCLUDE AN L.S TERM IS TAKEN, THE STARTING CONDITIONS FOR
C       THE INTEGRATION ARE SET UP AS IN HARTREE% CALCULATION
C       OF ATOMIC STRUCTURE, 81, OR FROM THE VALUE OF THE L.S TERM AT
C       THE SECOND MESH POINT BELOW WHICH A CONSTANT VALUE FOR THIS TERM
C       IS TAKEN TO AVOID DIVERGENCES.
      JR=0
      XV=DX/SPLIT
      IF(IREL.EQ.1)GOTO 65
      VO=-(ATP(1)-ATP(2))/DX
      GOTO 66
  65  VBCON=ALDS*(ATP(3)*X(2)-ATP(2)*X(3))/(X(2)**3*
     &      (X(3)-X(2)))/CC
      VO=VBCON*0.5
  66  ALPH=(2.0*Z*Z-(CL+1.0)*(2.0*VO+2.0*E))/
     &      (2.0*(2.0*CL+3.0)*(CL+1.0))
      BETA=Z*(2.0*Z*Z-(3.0*CL+4.0)*(2.0*VO+
     &     2.0*E))/(6.0*(CL+1.0)*(CL+2.0)*(2.0*CL+3.0))
      RL=1.0
      KDL2=L+1
      DO 2 K=1,KDL2
   2  RL=XV*RL
      RAD(1)=RL*(1.0-XV*Z/(CL+1.0)+ALPH*XV*XV-BETA*XV*XV*XV)
      RAD(2)=RL*(-Z/(CL+1.0)+2.0*ALPH*XV-3.0*beta*XV*XV)
      RAD(3)=XV
      N1=2
      INIT=0
      IF (M-1) 5,6,5
   5  IF (INIT) 3,10,3
   3  RAD(1)=ST1
      RAD(2)=ST2
      RAD(3)=ST3
  10  PRAD=STORE(JR+11)
   6  DO 11 K=1,10
      J2=IFIX(SPLIT)
      IF (J2-N1) 21,20,20
COMMENT INTEGRATION PROCEEDS OVER TEN VALUES OF THE BASIC GRID.
C       AFTER EACH RANGE OF TEN VALUES, AND BEFORE PROCEEDING
C       TO THE NEXT, THE STEP SIZE IS HALVED AND THE ACCURACY OF
C       THE CURRENT RANGE IS CHECKED. IF THE REQUIRED ACCURACY HAS
C       NOT BEEN REACHED AFTER HALVING THE STEP LENGTH 4 TIMES,
C       THE PROGRAM PRINTS OUT THE ACCURACY ACHIEVED, AND PROCEEDS.
C       THIS MONITOR OFTEN OCCURS IN THE FIRST RANGE OF TEN, AND
C       CAN BE IGNORED IF TEST IS SMALL. SUBROUTINE DA01A
C       ADVANCES THE INTEGRATION ACROSS ONE GRID INTERVAL
  20  DO 7 J=N1,J2
      XT=RAD(3)
      DDX=DX/SPLIT
      CALL DA01A(RAD,A,B,C,D,F,2,XT,DDX)
      RAD(3)=XT
      IF (BIG-ABS(RAD(1))) 12,12,7
  12  BIG=ABS(RAD(1))
   7  CONTINUE
  21  N1=1
  11  STORE(JR+K+1)=RAD(1)
      SPLIT=2.0*SPLIT
      M=M+1
      IF (M-5) 28,28,26
  26  TEST=ABS((PRAD-RAD(1))/BIG)
      WRITE (6,27) TEST
  27  FORMAT(' MONITOR M>5 TEST=',1pE13.5,'   Is it small?')
      GO TO 16
  28  IF (M-2) 13,14,13
  14  IF (INIT) 15,4,15
  15  IF (INIT-1) 13,5,13
  13  IF (ABS((PRAD-RAD(1))/BIG)-1.0E-4) 16,16,17
  17  IF (INIT) 18,4,18
  18  IF (INIT-1) 16,5,16
  16  ST1=RAD(1)
      ST2=RAD(2)
      ST3=RAD(3)
      SPLIT=1.0
      INIT=1
      IF (NF-JR-12) 23,23,22
  22  DX=2.0*DX
  23  JR=JR+10
      M=1
COMMENT THE PHASE SHIFT IS NOW CALCULATED AND SET IN DEL. DL IS
C       A LOGARITHMIC DERIVATIVE
      IF (RHO-X(JR+1)) 19,19,6
  19  DL=(RAD(2)+(CL+1.0)*RAD(1)/RAD(3))/RAD(1)-1.0/RAD(3)
      ES=SQRT(2.0*E)
      BFARG=RAD(3)*ES
      CALL CALCBF(BFARG,CL,BFJ1,BFN1)
      DCL=CL+1.0
      CALL CALCBF(BFARG,DCL,BFJ2,BFN2)
      DJL=CL*BFJ1/BFARG-BFJ2
      DNL=CL*BFN1/BFARG-BFN2
      A1=ES*DJL-DL*BFJ1
      A2=ES*DNL-DL*BFN1
      DEL=PI/2.0
      IF (ABS(A2)-1.0E-8) 25,25,24
  24  TD=A1/A2
      DEL=ATAN(TD)
C   CALCULATES ORIGIN OF PHASE SHIFT W.R.T. ZERO POTENTIAL
C   LOGARITHMIC DERIVATIVE AT NEAREST POINT OUTSIDE THE
C   MUFFIN-TIN RADIUS
C   FIRST CALCULATES NUMBER OF ZEROES IN WAVEFUNCTION
  25  NR=0
      DO 30 J=3,NMESH
      IF(STORE(J).LE.0.0.AND.STORE(J-1).GT.0.0)NR=NR+1
      IF(STORE(J).GE.0.0.AND.STORE(J-1).LT.0.0)NR=NR+1
      IF(J.EQ.JR+1)GOTO 41
  30  CONTINUE
  41  IF(L.GT.0)GOTO 35
      NZ=BFARG/PI
      GOTO 40
  35  AMU=1.0+4.0*CL*(CL+1.0)
      NZ=1
      AM1=AMU-1.0
C   CALCULATES APPROX. NO. OF ZEROES IN 'ZERO POTENTIAL' WAVEFUNCTION
C   BY A SERIES APPROXIMATION
  36  BT=(FLOAT(NZ)+CL/2.0)*PI
      BB=8.0*BT
      BZR=BT-AM1*(1.0/BB+4.0*(7.0*AMU-31.0)/3.0/BB**3
     &    +32.0*(83.0*AMU*AMU-982.0*AMU+3779.0)/15.0/BB**5
     &    +64.0*(6949.0*AMU*AMU*AMU-153855.0*AMU*AMU+1585743.0*AMU
     &    -6277237.0)/105.0/BB**7)
      IF(BFARG-BZR)38,39,37
C  CALCULATES NO. OF ZEROES IN DERIV. OF 'ZERO POTENTIAL' WAVEFTN.
  37  NZ=NZ+1
      GOTO 36
  38  NZ=NZ-1
C  CALCULATES EXACT NO. OF ZEROES IN 'ZERO POTENTIAL' WAVEFUNCTION
  39  IF(DJL.EQ.0.0)GOTO 50
      EVT=FLOAT(NZ)/2.0-FLOAT(NZ/2)
      IF(EVT.GT.0.1.AND.DJL.GT.0.0)NZ=NZ+1
      IF(EVT.LT.0.1.AND.DJL.LT.0.0)NZ=NZ+1
      IF(DJL.GT.0.0.AND.BFJ1.LT.0.0)NZ=NZ-1
      IF(DJL.LT.0.0.AND.BFJ1.GT.0.0)NZ=NZ-1
C  CORRECTION FOR INITIAL MODULO PI CALCULATION OF DEL
  40  IF(BFJ1.LT.0.0)GOTO 45
      IF(DEL.GT.0.0.AND.A1.LT.0.0)DEL=DEL-PI
      IF(DEL.LT.0.0.AND.A1.GT.0.0)DEL=DEL+PI
      GOTO 50
  45  IF(DEL.LT.0.0.AND.A1.LT.0.0)DEL=DEL+PI
      IF(DEL.GT.0.0.AND.A1.GT.0.0)DEL=DEL-PI
  50  DEL=DEL+FLOAT(NR-NZ)*PI
      ANORM=RAD(3)*(COS(DEL)*BFJ1-SIN(DEL)*BFN1)/RAD(1)
      DO 55 J=1,NMESH
  55  STORE(J)=STORE(J)*ANORM
      RETURN
      END




COMMENT SPHERICAL BESSEL FUNCTIONS 'J' AND 'N'
C       PENDRY CAVENDISH            JUNE 1970
      SUBROUTINE CALCBF(Z,CL,BJ,BN)
      L=IFIX(CL)
      IF (CL-0.5) 10,10,4
  10  BJ=1.0
      IF (1.0E-6-ABS(Z)) 11,11,8
  11  BJ=SIN(Z)/Z
      GO TO 8
   4  IF (1.1*CL-ABS(Z)) 12,12,6
  12  BJ=0.0
      IF (1.0E-6-ABS(Z)) 13,13,8
  13  BF2=SIN(Z)/Z
      BF1=SIN(Z)/(Z*Z)-COS(Z)/Z
      BFJ=2.0
      IF (L-1) 23,23,22
  22  DO 14 J=2,L
      BFS=BF1
      BF1=(2.0*BFJ-1.0)*BF1/Z-BF2
      BF2=BFS
  14  BFJ=BFJ+1.0
  23  BJ=BF1
      GO TO 8
   6  BFP=1.0
      BF=0.0
      BFLL=2.0*CL+3.0
      BFSR=1.0
      BFSX=0.5*Z*Z
   7  BF=BF+BFP
      BFP=-BFP*BFSX/(BFSR*BFLL)
      BFSR=BFSR+1.0
      BFLL=BFLL+2.0
      IF (ABS(BFP)-ABS(BF)*1.0E-6) 15,15,7
  15  BFLX=1.0
      BFNM=1.0
      DO 16 J=1,L
      BFLX=Z*BFLX
  16  BFNM=BFNM+1.0
      BFZ=1.0
      BFY=2.0*CL+1.0
      LL=IFIX(BFY)+1
      DO 17 JBB=1,LL,2
      BFZ=BFZ*BFY
  17  BFY=BFY-2.0
      BJ=BFLX*BF/BFZ
   8  IF (CL-0.5) 18,18,9
  18  BN=-COS(Z)/Z
      GO TO 5
   9  BF2=-COS(Z)/Z
      BF1=-COS(Z)/(Z*Z)-SIN(Z)/Z
      BFJ=2.0
      IF (L-1) 21,21,20
  20  DO 19 J=2,L
      BFS=BF1
      BF1=(2.0*BFJ-1.0)*BF1/Z-BF2
      BF2=BFS
  19  BFJ=BFJ+1.0
  21  BN=BF1
   5  RETURN
      END




COMMENT RUNGE-KUTTA-MERSON WITH ESTIMATED TRUNCATION ERROR E(I).
C       ADVANCES THE INTEGRATION OF THE SET OF DIFFERENTIAL
C       EQUATIONS  DY(I)/DX=F(I),I=1,N       BY ONE STEP FROM X0
C       TO X0+DX. THERE ARE FIVE ENTRIES TO DYBDX TO SET UP
C       F(J)=FUNC(Y(I),X),I=1,N
C       MCVICAR HARWELL                       FEBRUARY 1963
      SUBROUTINE DA01A(Y,E,A,B,C,D,N0,X0,DX)
      DIMENSION Y(3),E(2),A(2),B(2),C(2),D(2)
      N=N0
      Z=X0
      H=DX/3.
      DO 1 I=1,N
      D(I)=Y(I)
   1  CONTINUE
      X=Z
      DO 41 J=1,5
      CALL DYBDX(Y,E,N,X)
      DO 21 I=1,N
      GO TO (11,12,13,14,15),J
  11  A(I)=H*E(I)
      Y(I)=D(I)+A(I)
      GO TO 21
  12  B(I)=H*E(I)
      Y(I)=D(I)+(A(I)+B(I))*0.5
      GO TO 21
  13  B(I)=H*E(I)
      Y(I)=D(I)+(A(I)+B(I)*3.)*0.375
      GO TO 21
  14  C(I)=H*E(I)
      Y(I)=D(I)+(A(I)-B(I)*3.+C(I)*4.)*1.5
      GO TO 21
  15  Z=D(I)
      D(I)=H*E(I)
      Y(I)=Z+(A(I)+C(I)*4.+D(I))*0.5
      E(I)=(A(I)*2.-B(I)*9.+C(I)*8.-D(I))/10.
  21  CONTINUE
      GO TO (31,41,33,34,41),J
  31  X=Z+H
      GO TO 41
  33  X=Z+DX*0.5
      GO TO 41
  34  X=Z+DX
  41  CONTINUE
      X0=X
      RETURN
      END




COMMENT SUBROUTINE DYBDX SUPPORTS THE INTEGRATION SUBROUTINE DA01A
C       L.S TERM INCLUDED IF IREL=1.
      SUBROUTINE DYBDX(YT,F,N,XT)
      DIMENSION YT(2),F(2)
      COMMON /WS/K,JR,X(200),CL,E,ATP(200)
      COMMON /REL/IREL,ALDS,VBCON
      DX=X(JR+K+1)-X(JR+K)
      A=(XT-X(JR+K))/DX
      V=(ATP(JR+K)+(ATP(JR+K+1)-ATP(JR+K))*A)/XT
      IF(IREL.EQ.1)GOTO 5
      VBAR=0.0
      GOTO 10
   5  IF(JR+K.GE.2)GOTO 8
      VBAR=VBCON
      GOTO 10
   8  VBAR=ALDS*((ATP(JR+K+1)-ATP(JR+K))/DX-V)/18779.688/(XT*XT)
  10  F(1)=(YT(2)+(CL+1.0)*YT(1)/XT)
      F(2)=((2.0*V+VBAR-2.0*E)*YT(1)
     &     -(CL+1.0)*YT(2)/XT)
      RETURN
      END




COMMENT RONUC CALCULATES MATRIX ELEMENTS AND CROSS SECTIONS FOR ATOMIC EXAFS.
C       CALCULATES BOTH LI+1 AND LI-1 TRANSITIONS FOR ARBITRARY EDGES.
C       LI=INITIAL STATE ANGULAR MOMENTUM
C       L =FINAL STATE ANGULAR MOMENTUM
C       ncore = index(position) of initial state wf in the input file
C       ECORE=INITIAL STATE ENERGY : read in referenced to vacuum
C       CWM=INITIAL STATE WAVEFUNCTION
C       STORE=FINAL STATE WAVEFUNCTION
C       ANG=ANGULAR PART OF MATRIX ELEMENT
C       EMM=RADIAL PART OF MATRIX ELEMENT
C       PRM=EXAFS CROSS SECTION
      SUBROUTINE RONUC(Z,NM,VM,XM,WM,RM,E1,DE,NJE,LI,IR,PI)
      DIMENSION XM(260),VM(260),CWM(260),CWL(200)
      DIMENSION STORE(200),DING(200),PVPRRPHI(200)
      DIMENSION PRE(200),PREMM(200),PRPHS(200)
      COMMON /WS/KCF,JRCF,Y(200),CL,E,VL(200)
      COMMON /WF/WFC(260,25,4),WC(25),LC(25)

      READ(5,210) NCORE,ECORE
      IF(NCORE.LE.0)GOTO 31
C     Adjust core energy to be relative the muffin-tin zero level
      ECORE=ECORE-WM
      WRITE(6,211)NCORE,ECORE,LI

C     SET UP H-F GRID (HERMAN-SKILLMAN FORM).
c      110 pt mesh.  DX & Y = r values
c      NF should = 110 for nfinal of HS data, NMESH = 200 for extended mesh
      NMESH=200
      E=10.0
      AMU=0.88534138/EXP(ALOG(Z)/3.0)
      DX=0.01*AMU
      ES=SQRT(2.0*E)
      Y(1)=0.0
      J=2
      NF=NMESH
   21 KDL2=J+9
      DO 22 K=J,KDL2
   22 Y(K)=Y(K-1)+DX
      J=J+10
      IF(J-NMESH+10) 23,23,9
   23 IF(0.3/ES-DX) 5,5,24
   24 DX=DX+DX
      GO TO 21
    5 NF=J-1
      DO 7 K=J,NMESH
    7 Y(K)=Y(K-1)+DX
    9 CONTINUE

c     Read in the normalized initial state wavefunction.
c       Make the norm potential and norm wf flat beyond the muffin-tin radius.
c       Then, adjust the grid from Loucke's to HS.  Reset NMESH after each call
c       as CHGRID changes the value.
c       Finally, on HS grid, set norm pot and norm wf to zero beyond MT radius.
      CWM(1)=0.
      DO 1 J=2,NM
   1  CWM(J)=WFC(J-1,NCORE,IR)
      DO 13 J=1,NM
      IF (XM(J).GT.RM) CWM(J)=CWM(J-1)
  13  IF (XM(J).GT.RM) VM(J)=VM(J-1)
      CALL CHGRID(VM,XM,NM,VL,Y,NMESH)
      NMESH=200
      CALL CHGRID(CWM,XM,NM,CWL,Y,NMESH)
      NMESH=200
      DO J = 1,NMESH
        IF (Y(J).GT.RM) CWL(J) = 0.0
        IF (Y(J).GT.RM) VL(J)  = 0.0
      end do
c     Initialize values used in radial matrix element calculation that are
c       constant for both possible l_f
c       CWL = r*PHI for the initial state normalized wf
c       DV = d[VL]/dr  and VL = r*V
c        where d[VL] = V*dr + r*dV = V*dr + r*[(pV/pr)*dr] = [V + r*(pV/pr)]*dr
      DO J = 2,NMESH-2
        DV = 0.5*((VL(J)-VL(J-1))/(Y(J)-Y(J-1))+(VL(J+1)-VL(J))/
     &       (Y(J+1)-Y(J)))
        PVPRRPHI(J) = (DV-VL(J)/Y(J))/Y(J)*CWL(J)
      end do
      DING(1)=0.0
      WRITE(6,212)
      do J = 1, NMESH,2
        WRITE(6,213) Y(J),VL(J),CWL(J),Y(J+1),VL(J+1),CWL(J+1)
      end do


c     Start calculation of matrix element information.  
c     Begin with l_f = l_i + 1  and then repeat for l_f = l_i - 1
      L=LI+1

      DO 40 II=1,2
      WRITE (6,214) RM,Y(NMESH-3),L

c     CALCULATE ANGULAR part of matrix element
      XA = 0.0
      ZA = 0.0
      ANG = 0.0
      DO I = 1,2*L+1
        DO J = 1,2*LI+1
          MF = I-L-1
          MI = J-LI-1
          XA = BLM(L,-MF,LI,MI,1,1,PI)+BLM(L,-MF,LI,MI,1,-1,PI)
          ZA = BLM(L,-MF,LI,MI,1,0,PI)
          ANG = ANG+XA*XA+ZA*ZA
          WRITE(6,215)L,MF,LI,MI,XA,ZA
        end do
      end do
      ANG=4.0*PI/3.0*ANG
      FAC=64.0*PI*PI/9.0/137.03599976*ANG
      WRITE(6,216)ANG


c     Set up a loop over the energies requested to calculate the scattering
c     phase shifts and radial part of the matrix element for each energy
      E=E1
      WRITE(6,217)
      DO 18 JE=1,NJE
      CALL CALCFI(L,Z,STORE,DEL,NF,RM,PI)
c     Set up the integration variable DING =  (r*R) * [(pV/pr) * (r*PHI)]
c       STORE = r*R for the final state solution to schroed. eqtn
      DO 2 J=2,NMESH-2
        DING(J) = STORE(J)*PVPRRPHI(J)
   2  continue
c     Integrate out to the muffin-tin radius; Simpson's rule - step by 2
c       using natural progression along odds 1=init,3,5,7,... .  This will
c       take the integration past the break points where the step size changes
c       in a natural way.  110pt grid with 10 interval blocks of constant
c       step size; 1-11, 11-21, 21-31, ... .
      PREMM(JE) = 0.0
      DO 6 J = 3,NMESH-2,2
        if (Y(J) .lt. RM) then
          PREMM(JE) = PREMM(JE) +
     &                (DING(J-2) + 4.0*DING(J-1) + DING(J)) *
     &                (Y(J) - Y(J-1)) / 3.0
        endif
   6  continue
      PRE(JE) = E
      PRPHS(JE) = DEL
  18  E = E + DE

C     EMXS is the square of the radial matrix element divided by omega**2;
C       where omega is the angular frequency of the radiation causing the
C       transition (=E-ECORE) in Hartree atomic units. This quantity is also
C       written onto unit 8 in a format suitable for input to the XANES
C       program (this form of the matrix element is correctly normalised for
C       this purpose).
C     EMMM is the EXAFS cross section from the radial matrix element
      WRITE(6,218)
      do J = 1, NJE
        OMEGA = PRE(J) - ECORE
        EMXS=PREMM(J)**2/(OMEGA)**2
        EMMM=FAC*SQRT(2.0*PRE(J))*EMXS/(OMEGA)
        WRITE (6,219) PRE(J),PRPHS(J),OMEGA,PREMM(J),EMXS,EMMM
        WRITE(8,220) PRE(J),EMXS
      end do

      IF(LI.EQ.0)GOTO 31
      L=L-2
  40  CONTINUE

  31  RETURN
 210  FORMAT(I4,F11.5)
 211  FORMAT(/,'   INITIAL STATE  ::  Ncore=',I3,'   Ecore=',F11.5,
     &       '   LI=',I3)
 212  FORMAT(/,'   R*V(R), AU*HARTREES, AND R*PSI_CORE ON Her-Skil',
     &         ' 10pt-block GRID',/,
     &         1x,2('   R           R*V(R)      R*PSI(R) '))
 213  FORMAT(1x,6(1pE12.4))
 214  FORMAT(/,'   FINAL STATE  ::  MTRad=',F8.4,
     &         '   Rmax=',F8.4,'   L=',I3,//,
     &         '   ANGULAR MATRIX ELEMENTS',/,
     &         3x,'Lf',2x,'Mf',2x,'Li',2x,'Mi',5x,'XA',9x,'ZA')
 215  FORMAT(1x,4I4,2F11.5)
 216  FORMAT('              ANGULAR PART',F9.5)
 217  FORMAT(/,'   RADIAL MATRIX ELEMENTS AND PHASE SHIFTS ')
 218  FORMAT(5X,'E',8X,'PHS',6X,'OMEGA',5X,'EMM',8x,'(EMM/omega)^2',
     &       2X,'Cross Section')
 219  FORMAT(1x,2(F9.4),f10.4,1pE13.5,1pe13.5,1pe15.5)
 220  FORMAT(2(1pE12.5))
      END




COMMENT FUNCTION BLM PROVIDES THE INTEGRAL OF THE PRODUCT
C       OF THREE SPHERICAL HARMONICS, EACH OF WHICH CAN BE
C       EXPRESSED AS A PRE FACTOR TIMES A LEGENDRE FUNCTION.
C       THE INTEGRAL OF THE THREE LEGENDRE FUNCTIONS FOLLOWS
C       GAUNT'S SUMMATION SCHEME SET OUT BY SLATER% ATOMIC
C       STRUCTURE, VOL1, 309,310
      FUNCTION BLM(IL1,IM1,IL2,IM2,IL3,IM3,PI)

      BLM=0.0
      L1=IL1
      L2=IL2
      L3=IL3
      M1=IABS(IM1)
      M2=IABS(IM2)
      M3=IABS(IM3)
      IF(L1.GE.M1.AND.L2.GE.M2.AND.L3.GE.M3)GOTO 2
      WRITE(6,230)IL1,IM1,IL2,IM2,IL3,IM3
      RETURN
   2  M=IABS(IM1+IM2+IM3)+MOD(L1+L2+L3,2)
      IF(M.NE.0.OR.L1.GT.(L2+L3).OR.L1.LT.IABS(L2-L3))RETURN
      L=L1
      M=M1
      M1=MAX0(M1,M2,M3)
      IF(M1.EQ.M)GOTO 5
      IF(M1.EQ.M2)GOTO 3
      IF(M1.EQ.M3)GOTO 4
   3  L1=L2
      L2=L
      M2=M
      GOTO 5
   4  L1=L3
      L3=L
      M3=M
   5  IF(L2.GE.L3)GOTO 6
      L=L2
      M=M2
      L2=L3
      M2=M3
      L3=L
      M3=M
   6  IS=(L1+L2+L3)/2
      BLM=(-1.0)**(IS-L2-M3+M1)*SQRT(FLOAT((2*L1+1)*(2*L2+1)*(2*L3+1))
     &   /PI)/FLOAT(2*(2*IS+1))*SQRT(FAC(L1+M1,M1+M1))/FAC(IS-L1,IS-L1)
     &   *SQRT(FAC(L2+M2,M2+M2))/FAC(IS-L2,IS-L2)*SQRT(FAC(L3+M3,M3+M3))
     &   *FAC(MAX0(L2+L3-M1,L2-L3+M1),IABS(2*(L3-M1)))**ISIGN(1,L3-M1)
      IF(L3.EQ.0)RETURN
      DO 7 I=1,L3
   7  BLM=BLM/FLOAT(2*(2*(IS-I)+1))
      A=1.0
      AA=1.0
      I1=L1+M1
      I2=MAX0(1,L2+L3-M1)
      I3=0
      I4=L1-M1
      I5=L2-L3+M1
      I6=L3-M3
      IT=MIN0(I4,I6)
      IF(IT.EQ.0)RETURN
      DO 8 I=1,IT
      I1=I1+1
      I3=I3+1
      I5=I5+1
      AA=-AA*FLOAT(I1*I4*I6)/FLOAT(I2*I3*I5)
      A=A+AA
      I2=I2-1
      I4=I4-1
   8  I6=I6-1
      BLM=A*BLM
      RETURN
 230  FORMAT(' *** INVALID ARGUMENTS FOR BLM ',5(I2,','),I2)
      END




      FUNCTION FAC(N,M)
      FAC=1.0
      IF(M.EQ.0)RETURN
      M1=N-M+1
      DO 1 I=M1,N
   1  FAC=FAC*FLOAT(I)
      RETURN
      END




      SUBROUTINE CALPFI(FI,DFI,SI,DSI,DEL,EE,L,ISPL,INDP,IP,X,ATP,NMESH,
     &                  KRMT,Z,NF)
C     =================
C
C     RADIAL INTEGRATION BY A PREDICTOR-CORRECTOR METHOD.
C     SPLITS S.E. INTO 4 REAL COUPLED FIRST ORDER ODE'S
C
C     INPUTS
C     ======
C
C     E         ENERGY
C     L         ANGULAR MOMENTUM
C     ISPL      INTEGRATION PARAMETER (1.0)
C     INDP      IF = 0 WILL CALCULATE REGULAR SOLN. & PHASE SHIFTS
C               IF =-1 WILL CALCULATE IRREGULAR SOLUTION
C     NMESH     # OF MESH POINTS ON H/S GRID
C     KRMT      MESH POINT BEYOND MUFFIN-TIN
C     NF        MESH POINT BEYOND WHICH GRID DOUBLING CEASES
C     Z         ATOMIC NUMBER
C     X(I)      ITH H/S GRID POINT
C     ATL(I)    -R*V(R) AT ITH H/S GRID POINT
C
C     OUTPUTS
C     =======
C
C     FI(I)     REGULAR WAVEFUNCTION AT ITH H/S GRID POINT
C     DFI(I)    DERIVATIVE OF REGULAR WAVEFUNCTION AT ITH H/S...
C     SI(I)     IRREGULAR WAVEFUNCTION AT ITH H/S GRID POINT
C     DSI(I)    DERIVATIVE OF IRREGULAR WAVEFUNCTION AT ITH H/S...
C     DEL       PHASE SHIFT
C
C     SETTING INDP=-1 PRODUCES THE IRRGULAR WAVEFUNCTION IN SI & DSI
C     INTEGRATION PROCEDES FROM MUFFIN TIN RADIUS INWARDS
C     IF INDP.NG.0 THE SUBROUTINE MAY BE CALLED WITH
C     ARGUMENTS -           (FI,DFI,FI,DFI,....  0,IP)
C     METHOD USES SET OF COUPLED FIRST ORDER EQUATIONS DERIVED
C     FROM THE SECOND ORDER SCHROEDINGER EQUATION@D
C     R1=FI                     R2=DFI-(L+1)/R *FI
C     DERIVATIVES ARE CALCULATED IN DYBDX@D
C       R1'=(L+1)/R *R1 +R2      R2'=-(2*(E-V))*R1 -(L+1)/R *R2
C
      COMPLEX FI,DFI,SI,DSI,DEL,E,EE,RN
      COMPLEX ST1,ST2,CZ,ZZ,CI
      COMPLEX BFJ1,BFJ2,BFN1,BFN2,H1,H2,DH1,DH2
      COMPLEX ALPH,BETA,GAMMA,A1,A2,A3,VO
      DIMENSION FI(200),DFI(200),SI(200),DSI(200)
      DIMENSION X(200),ATP(200),RN(42),ND(7)
      DATA NR/10/,  CZ,CI/(0.0,0.0),(0.0,1.0)/
      DATA  ND/1,7,13,19,25,31,37/
C
      ABSQ(ZZ)=REAL(ZZ)**2+AIMAG(ZZ)**2
C
C     INITIALISE FOR REGULAR SOLUTION
C
      E=EE
   10 DO 100 J=1,NMESH
      SI(J)=CZ
 100  DSI(J)=CZ
C
      CL=FLOAT(L)
      IREG=1
      SPL=ISPL
      DX=X(2)-X(1)
C
      VO=E+ (3.0*ATP(1)-4.0*ATP(2)+ATP(3))/(2.0*DX)
      V1=-(ATP(1)-2.0*ATP(2)+ATP(3))/(2.0*DX*DX)
      A1=-Z/(CL+1.0)
      A2=(Z*Z-(CL+1.0)*VO)/((2.0*CL+3.0)*(CL+1.0))
      A3=-(Z*(Z*Z-(3.0*CL+4.0)*VO)+(CL+1.0)*(2.0*CL+3.0)*V1)/
     &    (3.0*(CL+1.0)*(CL+2.0)*(2.0*CL+3.0))
C
C     SET THE POWER SERIES EXPANSION TO THIRD ORDER IN THE FIRST THREE
C     POINTS FROM THE ORIGIN AS STARTING VALUES FOR THE PREDICTOR
C     -CORRECTOR, THESE ARE@D-
C
C          FI=R**(L+1) *(A0 +A1*R +A2*R**2 +A3*R**3 ...)
C          R1=R**(L+1) *(A0 +A1*R +A2*R**2 +A3*R**3 ...)
C          R2=R**(L+1) *(A1 +2*A2*R +3*A3*R**2 ...)
C
      K=ND(4)-1
      DO 110 J=1,4
 110  RN(K+J)=CZ
C
      IF(L.LE.0) THEN
       RN(K+3)=1.0
       RN(K+4)=A1
      ENDIF
C
      JRK=1
      JR=NR
      R=0.0
      DX=DX/SPL
      JSPL=ISPL
C
      DO 16 J=2,4
      JSPL=JSPL-1
      R=R+DX
      B=R
      K=0
C
   14 IF (K.NE.L) THEN
         K=K+1
         B=R*B
         GOTO 14
      ENDIF
C
      K=ND(5-J)
      RN(K)=B*(1.0+R*(A1+R*(A2+R*A3)))
      RN(K+1)=B*(A1+R*(2.0*A2+R*(3.0*A3)))
      CALL DYBDX1(X,ATP,CL,EE,NMESH,KRMT,JRK,RN(K),RN(K+2),2,R)
      RN(K+4)=CZ
      RN(K+5)=CZ
C
      IF (JSPL.LE.0) THEN
         JSPL=ISPL
         JR=JR-1
         JRK=JRK+1
         SI(JRK)=RN(K)
         DSI(JRK)=RN(K+1)+(CL+1.0)*RN(K)/R
      ENDIF
C
 16   CONTINUE
C
C     INTEGRATION PROCEEDS OUT TO RMT
C     IT IS ASSUMED THAT THE POINTS ARE ON HERMAN-SKILLMAN GRID
C
 20   CONTINUE
C
      DO 23 J=1,JR
 21   JSPL=JSPL-1
      R=R+DX
      NK=ND(7)
      DO 22 K=2,7
 22   ND(9-K)=ND(8-K)
      ND(1)=NK
C
      CALL DA02A(RN(ND(1)),RN(ND(2)),RN(ND(3)),RN(ND(4)),RN(ND(5)),
     &           4,R,DX,X,ATP,CL,EE,NMESH,KRMT,JRK)
C
      IF (JSPL.GT.0) GO TO 21
      JSPL=ISPL
      IF(IREG.GT.0) JRK=JRK+IREG
      SI(JRK)=RN(ND(1))
      DSI(JRK)=RN(ND(1)+1)+(CL+1.0)*RN(ND(1))/R
      R=X(JRK)
      IF(IREG.LT.0) JRK=JRK+IREG
 23   CONTINUE
C
      IF(JRK.GE.KRMT .OR. JRK.LE.1) GO TO 30
      DX=(X(JRK+IREG)-X(JRK))/SPL
      JR=NR
      IF(JRK.GE.NF) GO TO 20
      IF(IREG.LT.0) GO TO 24
      NK=ND(2)
      ND(2)=ND(3)
      ND(3)=ND(5)
      ND(5)=ND(4)
      ND(4)=ND(7)
      ND(7)=NK
      GO TO 20
C
 24   DO 25 J=1,4
      K=J-1
      RN(ND(6)+K)=(35*RN(ND(1)+K)+140*RN(ND(2)+K)-70*RN(ND(3)+K)
     &            +28*RN(ND(4)+K)-5*RN(ND(5)+K))/128.0
      RN(ND(7)+K)=(-5*RN(ND(1)+K)+60*RN(ND(2)+K)+90*RN(ND(3)+K)
     &            -20*RN(ND(4)+K)+3*RN(ND(5)+K))/128.0
 25   CONTINUE
C
      NK=ND(2)
      ND(2)=ND(6)
      ND(6)=ND(4)
      ND(4)=ND(7)
      ND(7)=ND(5)
      ND(5)=ND(3)
      ND(3)=NK
      IF(JRK.EQ.NR) JR=NR-1
      GO TO 20
C
C     THE PHASE SHIFT IS NOW CALCULATED AND SET IN DEL. DL IS
C     A LOGARITHMIC DERIVATIVE.  BOTH ARE EVALUATED AT THE GRID
C     POINT IMMEDIATELY BEYOND RMT. FI HOLDS VALUES OF THE REGULAR
C     WAVE FUNCTION TIMES RADIUS, AND DFI HOLDS DERIVATIVES OF FI
C
 30   IF(IREG.LT.0) GO TO 50
      JMAX=KRMT-1
 31   B=ABSQ(SI(JMAX))+ABSQ(DSI(JMAX))
      JMAX=JMAX+1
      IF(ABSQ(SI(JMAX)).LT.B*1.E-8 .AND. JMAX.LT.JRK) GOTO 31
      ST1=DSI(JMAX)/SI(JMAX)-1.0/X(JMAX)
      ST2=SI(JMAX)/X(JMAX)
      BETA=CSQRT(2.0*E)
      ALPH=X(JMAX)*BETA
      IF(ABS(AIMAG(ALPH)).GT.5.0) GO TO 34
      CALL CALXBF(ALPH,CL,BFJ1,BFN1)
      CALL CALXBF(ALPH,CL+1.0,BFJ2,BFN2)
      BFJ2=BETA*(CL*BFJ1/ALPH-BFJ2)
      BFN2=BETA*(CL*BFN1/ALPH-BFN2)
      H1=BFJ1+CI*BFN1
      DH1=BFJ2+CI*BFN2
      GAMMA=(BFJ2-ST1*BFJ1)/(BFN2-ST1*BFN1)
      IF(ABSQ(GAMMA).GT.1.E-2) GO TO 32
      ALPH=GAMMA*GAMMA
      DEL=GAMMA*(1.0-ALPH*(1.0/3.0-ALPH*(1.0/5.0-ALPH*(1.0/7.0))))
      ID=1
      GO TO 33
 32   ALPH=CI*GAMMA
      DEL=CMPLX(0.0,-0.5)*CLOG((1.0+ALPH)/(1.0-ALPH))
      ID=2
 33   ALPH=CCOS(DEL)*(BFJ1-GAMMA*BFN1)/ST2
      GAMMA=CEXP(CI*DEL)
      GO TO 35
 34   CALL CALXBH(ALPH,CL,H1,H2)
      CALL CALXBH(ALPH,CL+1.0,DH1,DH2)
      DH1=BETA*(CL*H1/ALPH-DH1)
      DH2=BETA*(CL*H2/ALPH-DH2)
      DEL=CMPLX(0.0,-0.5)*CLOG((DH2-ST1*H2)/(ST1*H1-DH1))
      ID=3
      GAMMA=CEXP(CI*DEL)
      ALPH=0.5*(H1*GAMMA+H2/GAMMA)/ST2
C
C     FI AND DFI ARE NORMALISED AT THE MUFFIN-TIN RADIUS
C
 35   DO 36 J=1,NMESH
      FI(J)=ALPH*SI(J)
 36   DFI(J)=ALPH*DSI(J)
      IF(INDP.GT.0) GO TO 40
      BETA=CZ
      GAMMA=CZ
      GO TO 60
C
C       INITIALISE FOR IRREGULAR SOLUTION
C
 40   IREG=-1
      H1=H1*GAMMA
      DH1=DH1*GAMMA
      H2=ALPH*ST2
      DH2=H2*ST1
C
      DO 41 J=1,4
      K=ND(J)
      R=R+DX
      ST2=R*BETA
      CALL CALXBH(ST2,CL,BFJ1,BFN1)
      CALL CALXBH(ST2,CL+1.0,BFJ2,BFN2)
      ST1=BFJ1*GAMMA
      ST2=BETA*(CL*BFJ1/ST2-BFJ2)*GAMMA
      RN(K)=ST1
      RN(K+1)=ST2-(CL+1.0)*ST1/R
      CALL DYBDX1(X,ATP,CL,EE,NMESH,KRMT,JRK,RN(K),RN(K+2),2,R)
      RN(K+4)=CZ
      RN(K+5)=CZ
 41   CONTINUE
C
      R=X(JRK)+DX
      DX=-DX
      JSPL=0
      JR=NR+1
      GO TO 20
C
C       NORMALISE IRREGULAR SOLUTION TO RETURN FI+  IN SI,DSI
C
 50   ST1=DSI(JMAX)/SI(JMAX)-1.0/X(JMAX)
      ST2=SI(JMAX)/X(JMAX)
      BETA=(H1*DH2-DH1*H2)/((DH2-ST1*H2) *ST2)
      GAMMA=(DH1-ST1*H1)/(DH2-ST1*H2)
C
      DO 51 J=1,NMESH
      SI(J)=BETA*SI(J) +GAMMA*FI(J)
 51   DSI(J)=BETA*DSI(J) +GAMMA*DFI(J)
C
C       WRITE REQUIRED OUTPUT
C
c60   CL=RMT
 60   CL=KRMT
      IF(IP.LE.3) GO TO 64
      JR=(JMAX+1)/2
      WRITE(6,240)E,L,ID
      WRITE(6,241) ALPH,BETA,GAMMA,X(JMAX)
      WRITE(6,242) ((X(K),FI(K),DFI(K),K=J,JMAX,JR),J=1,JR)
      IF(INDP.EQ.0) GO TO 64
      WRITE(6,242) ((X(K),SI(K),DSI(K),K=J,JMAX,JR),J=1,JR)
   64 RETURN
  240 FORMAT(' ENERGY = ',F7.3,5X,'L =',I3,'  ID =',I3)
  241 FORMAT(' ALPH,BETA,GAMMA =',3(E14.4,E12.4),5X,'AT RADIUS',E12.4)
  242 FORMAT(' X(K),FI(K),DFI(K)',/,2(E15.5,2(E13.5,E12.5)))
      END




      SUBROUTINE DA02A(A,B,C,D,E,N,X,DX,XR,ATP,CL,EE,NMESH,KRMT,JRK)
C     ================
      COMPLEX XR(1),ATP(1),EE
C
C     SUBROUTINE DA02A IS A PREDICTOR-CORRECTOR INTEGRATION ROUTINE
C     FOR A SET OF COUPLED FIRST ORDER DIFERRENTIAL EQUATIONS.
C     SEE HAMMING, CH 15.
C     A RETURNS  THE FUNCTIONS AND DERIVATIVES AT X
C     B,C...  CONTAIN THE FUNCTIONS AND DERIVATIVES AT X-DX,X-2DX...
C     EACH IS EFFECTIVELY DIMENSIONED  (N,3)%
C        (*,1)       FNS,  (*,2)  DERIVS,  (*,3)  (P-C).
C     DYBDX  SHOULD RETURN THE DERIVS  GIVEN THE FNS AT X.
C
      DIMENSION A(N),B(N),C(N),D(N),E(N)
      DATA A1,A2/0.0,0.125/,  N0/1/
      IF (N0.EQ.0) GO TO 1
      A0=1.0-A1-A2
      B0=55.0+9.0*A1+8.0*A2
      B1=-59.0+19.0*A1+32.0*A2
      B2=37.0-5.0*A1+8.0*A2
      B3=-9.0+A1
      DM=9.0-A1
      D0=19.0+13.0*A1+8.0*A2
      D1=-5.0+13.0*A1+32.0*A2
      D2=1.0-A1+8.0*A2
      AM1=(251.0-19.0*A1-8.0*A2)/(270.0-30.0*A1)
      AM2=(-19.0+11.0*A1-8.0*A2)/(270.0-30.0*A1)
      N0=0
    1 CONTINUE
      H=DX/24.0
      NN=N+N
      DO 2 J=1,N
      A(J+NN)=A0*B(J)+A1*C(J)+A2*D(J)
     &        +H*(B0*B(J+N)+B1*C(J+N)+B2*D(J+N)+B3*E(J+N))
    2 A(J)=A(J+NN)-AM1*C(J+NN)
      CALL DYBDX1(XR,ATP,CL,EE,NMESH,KRMT,JRK,A(1),A(N+1),N,X)
      DO 3 J=1,N
      CC=A0*B(J)+A1*C(J)+A2*D(J)
     &   +H*(DM*A(J+N)+D0*B(J+N)+D1*C(J+N)+D2*D(J+N))
      A(J+NN)=A(J+NN)-CC
    3 A(J)=CC+AM2*C(J+NN)
      CALL DYBDX1(XR,ATP,CL,EE,NMESH,KRMT,JRK,A(1),A(N+1),N,X)
      RETURN
      END




C     SUBROUTINE DYBDX SUPPORTS THE INTEGRATION SUBROUTINE DA01A
      SUBROUTINE DYBDX1(X,ATP,CL,EE,NMESH,KRMT,JRK,YT,F,N,XT)
C     ================
C
      COMPLEX EE,YT,F
      DIMENSION YT(2),F(2),ATP(1),X(1)
C
      V=0.0
C
C THE NEXT LINE WAS IN THE PROGRAM I INHERITED. IT DOES NOT MAKE SENSE.
C IT HAS BEEN REPLACED BY THE LINE FOLLOWING IT. THE NEW VERSION OF THE
C PROGRAM PRODUCES THE SAME OUTPUT AS THE OLD ONE, IMPLYING THAT JJ HAD
C BEEN INITIALIZED TO ZERO BY DEFAULT.
C  THE OLD LINE CONTAINED: J=JJ
C  THIS HAS BEEN REPLACED BY THE FOLLOWING LINE (ON 15 FEB. 1990):
      J=0
C
      IF(J.LE.0) J=1
      IF(J.GT.NMESH-2) J=NMESH-2
   10 IF(XT.GT.X(J) .OR. J.LE.1) GO TO 11
      J=J-1
      GO TO 10
   11 IF(J.GE.KRMT) GO TO 13
      IF(XT.LE.X(J+1) .OR. J.GE.NMESH-2) GO TO 12
      J=J+1
      GO TO 11
   12 X1=XT-X(J)
      X2=XT-X(J+1)
      X3=XT-X(J+2)
      D1=X2-X3
      D2=X3-X1
      D3=X1-X2
      V=-(D1*X2*X3*ATP(J)+D2*X3*X1*ATP(J+1)+D3*X1*X2*ATP(J+2))/
     &   (D1*D2*D3 *XT)
   13 F(1)=(CL+1.0)/XT *YT(1) +YT(2)
      F(2)=2.0*(V-EE) *YT(1) -(CL+1.0)/XT *YT(2)
      RETURN
      END




      SUBROUTINE CALXBF(Z,CL,BJ,BN)
C      ================
      COMPLEX ZZ, Z,BJ,BN, BF,BF1,BF2,BFP,BFS,BFSX,BFLX
      ABSQ(ZZ)=REAL(ZZ)**2+AIMAG(ZZ)**2
      IND=0
      L=IFIX(CL)
      IF (CL-0.5) 10,10,4
  10  BJ=1.0
      IF (1.0E-12-ABSQ(Z)) 11,11,8
  11  BJ=CSIN(Z)/Z
      GO TO 8
   4  IF ((1.1*CL)**2-ABSQ(Z)) 12,12,6
  12  BJ=0.0
      IF (1.0E-12-ABSQ(Z)) 13,13,8
  13  BF2=CSIN(Z)/Z
      BF1=CSIN(Z)/(Z*Z)-CCOS(Z)/Z
      BFJ=2.0
      IF (L-1) 23,23,22
  22  DO 14 J=2,L
      BFS=BF1
      BF1=(2.0*BFJ-1.0)*BF1/Z-BF2
      BF2=BFS
  14  BFJ=BFJ+1.0
  23  BJ=BF1
      GO TO 8
   6  BFP=1.0
      BF=0.0
      BFLL=2.0*CL+3.0
      BFSR=1.0
      BFSX=0.5*Z*Z
   7  BF=BF+BFP
      BFP=-BFP*BFSX/(BFSR*BFLL)
      BFSR=BFSR+1.0
      BFLL=BFLL+2.0
      IF (ABSQ(BFP)-ABSQ(BF)*1.0E-12) 15,15,7
  15  BFLX=1.0
      BFNM=1.0
      DO 16 J=1,L
      BFLX=Z*BFLX
  16  BFNM=BFNM+1.0
      BFZ=1.0
      BFY=2.0*CL+1.0
      LL=IFIX(BFY)+1
      DO 17 JBB=1,LL,2
      BFZ=BFZ*BFY
  17  BFY=BFY-2.0
      BJ=BFLX*BF/BFZ
   8  IF (CL-0.5) 18,18,9
  18  BN=-CCOS(Z)/Z
      GO TO 5
   9  BF2=-CCOS(Z)/Z
      BF1=-CCOS(Z)/(Z*Z)-CSIN(Z)/Z
  30  BFJ=2.0
      IF(L-1) 21,21,20
  20  DO 19 J=2,L
      BFS=BF1
      BF1=(2.0*BFJ-1.0)*BF1/Z-BF2
      BF2=BFS
  19  BFJ=BFJ+1.0
  21  BN=BF1
      IF(IND.GT.0) GO TO 31
   5  RETURN
      ENTRY        CALXBH(Z,CL,BJ,BN)
C     =================
C
      IND=1
      L=IFIX(CL)
      BF=CMPLX(0.0,1.0)*Z
 32   BF1=CEXP(BF)/BF
      IF(L.LE.0) GOTO 21
      BF2=BF1
      BF1=BF1*(1.0-BF)/Z
      IF(L.EQ.1) GOTO 21
      GOTO 30
 31   IND=0
      BJ=BN
      BF=-BF
      GOTO 32
      END



__________________________________________________________________________
Go to top of page