CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  program sumpro
C  written by Andreas Eichler (band)
C  modifyed by Georg Kern 
C
C  topic: finding surface-states
C
C  input: PROCAR
C         SUMPROCTR:
C                   1 2           ... ions for criterium 1
C                   .75 .50 .25   ... % for criterium 1
C                   1 2 6 8 4     ... ions for criterium 2
C                   .75 .50 .30   ... % for criterium 2
C                   5.3342        ... E(fermi) (optional)
C                   64            ... nr. of electrons (optional)
C
C  output: SBAND.1 ... for mathematica (sband.m) criterium 1 
C          SBAND.2 ... detto, but criterium 2
C          SBAND.user (for user - easy to read)
C
C  parameter: NKDIM ... maximum k-points
C             NBDIM ... maximum bands
C             NDIM  ... maximum ions
C
C  anything else is (nearly) selfexplaining
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PROGRAM SUMPRO
      PARAMETER(NKDIM=128,NBDIM=512,NDIM=256)
      CHARACTER*256 LINE
      REAL RELEN(NKDIM,NBDIM),PARSUM(NKDIM,NBDIM,NDIM)
      REAL CRI(2,3),OCCUP
      INTEGER ION(2,0:NDIM)
      LOGICAL LPROFM
      DIMENSION PSUM(2)	

      OPEN(10,FILE='PROCAR',STATUS='OLD')
      OPEN(11,FILE='SBAND.1')
      OPEN(12,FILE='SBAND.user')
      OPEN(13,FILE='SBAND.2')
      READ(10,'(A)') LINE
      IF (LINE(1:10).EQ.'PROCAR new') THEN
        LPROFM=.TRUE.
        READ(10,'(A)') LINE
      ELSE
        LPROFM=.FALSE.
      ENDIF

C here comes the first format for PROCAR
      READ(LINE(15:20),*,ERR=99) NKPTS
      READ(LINE(33:38),*,ERR=99) NBANDS
      READ(LINE(49:55),*,ERR=99) NIONS
      GOTO 100
C and now the second
 99   READ(LINE(15:19),*) NKPTS
      READ(LINE(40:43),*) NBANDS
      READ(LINE(63:66),*) NIONS

 100  WRITE(*,*) '   NIONS NKPTS NBANDS'
      WRITE(*,'(3I6)') NIONS,NKPTS,NBANDS
      IF ((NKPTS.GT.NKDIM).OR.(NBANDS.GT.NBDIM).OR.(NIONS.GT.NDIM)) THEN
        WRITE(*,*) 'ERROR: FIELDS TOO SMALL!!!'
        WRITE(*,*) '       CHECK NKDIM, NBDIM, NDIM'
        WRITE(*,*) 'PROGRAMM WILL BE STOPPED!!!!'
        GOTO 2401
      ENDIF

      OPEN(9,FILE='SUMPROCTR',STATUS='OLD',ERR=2303)
      GOTO 2310
 2303 PRINT*,'ERROR: SUMPROCTR DOES NOT EXIST!!!'
      PRINT*,'PROGRAMM WILL BE STOPPED!!!!'
      GOTO 2401
 2310 DO I=1,2
        READ(9,'(A)') LINE
        N=NROFI(LINE)
        IF (N.GT.NDIM) THEN 
          PRINT*,'TOO MANY IONS IN SUMPROCTR (cri ',I,')!!!'
          PRINT*,'PROGRAMM WILL BE STOPPED!!!!'
          CLOSE(9)
          GOTO 2401
        ENDIF
        READ(LINE,*) (ION(I,II),II=1,N)
        ION(I,0)=N
        READ(9,*) (CRI(I,II),II=1,3)
      ENDDO
      READ(9,*,END=2309,ERR=2309) EFERMI
      GOTO 2308
 2309 EFERMI=-1.e10
 2308 READ(9,*,END=2307,ERR=2307) NREL
      GOTO 2306
 2307 NREL=-1
 2306 DO I=1,2
        PRINT*,'criterium ',I
        PRINT*,'   ions: ',(ION(I,II),II=1,ION(I,0))
        PRINT*,'      %: ',CRI(I,:)
      ENDDO
      IF(EFERMI<-1.e9) THEN
        PRINT*,'E(fermi) will be calculated'
      ELSE
        PRINT*,'E(fermi) = ',EFERMI
      ENDIF
      IF(NREL>0) THEN
        PRINT*,'nr. of e- = ',NREL
      ELSE
        PRINT*,'nr. of e- not given'
        PRINT*,' (-> not very good calculated E(fermi))'
      ENDIF
      CLOSE(9)

      EFTEST=-1.e10
      WRITE(*,'(A,$)') 'reading PROCAR'
      DO 2300 NK=1,NKPTS
        WRITE(*,'(A,$)') '.'
        IF (LPROFM) THEN
         DO I=1,3
          READ(10,'(A)') LINE
         ENDDO
        ELSE
         DO I=1,4
          READ(10,'(A)') LINE
         ENDDO
        ENDIF

        DO 2300 NB=1,NBANDS
          NI=1
          IF (LPROFM) THEN
           READ(10,3203) NDUM,RELEN(NK,NB),OCCUP
 3203 FORMAT('band ',I3,' # energy',F14.8,' # occ.',F12.8//
     &       ' ion   s      p      d     tot')
          ELSE
           READ(10,*) NDUM,RELEN(NK,NB),OCCUP,NDUM,
     &      (RDUM,I=1,3),PARSUM(NK,NB,NI)
          ENDIF
C--this calculation method for the Fermi energy is not very good
C  only if NREL not given
          IF ((OCCUP.GT.1).AND.(RELEN(NK,NB).GT.EFTEST)) 
     &      EFTEST=RELEN(NK,NB)
C--don't know why next if block
C          IF (NIONS.G1) THEN

          IF (LPROFM) THEN
            DO NI=1,NIONS
              READ(10,*) NDUM,(RDUM,I=1,3),PARSUM(NK,NB,NI)
            ENDDO
            READ(10,'(A)') LINE
            READ(10,'(A)') LINE
          ELSE
 	    DO NI=2,NIONS
 	      READ(10,*) NDUM,(RDUM,I=1,3),PARSUM(NK,NB,NI)
            ENDDO
            READ(10,'(A)') LINE
            READ(10,'(A)') LINE
          ENDIF

C          ENDIF
2300    CONTINUE
        WRITE(*,*)

C--now calculate the Fermi energy with NREL
       IF (NREL>0) THEN
         NSOLL=NREL*NKPTS/2
         EFTEST2=MINVAL(RELEN)
         DELTAE=1
         DO WHILE (DELTAE>=0.00005)
           NIST=COUNT(RELEN(1:NKPTS,1:NBANDS).LE.EFTEST2)
           DO WHILE (NIST<NSOLL) 
             EFTEST2=EFTEST2+DELTAE
             NIST=COUNT(RELEN(1:NKPTS,1:NBANDS)<=EFTEST2)
           ENDDO
C           PRINT*,'DELTAE=',DELTAE,', EFERMI=',EFTEST2,', ist=',
C     &            NIST,', soll=',NSOLL
           EFTEST2=EFTEST2-DELTAE
           DELTAE=DELTAE/10
         ENDDO
       ELSE
         EFTEST2=EFTEST
         PRINT*,'Warning: calculated E(fermi) uncertain'
       ENDIF

C--subtract Fermi energy
       IF (EFERMI.LT.-1e9) THEN
         PRINT*,'using caclulatet fermi-energy'
         RELEN=RELEN-EFTEST2
       ELSE
         PRINT*,'using given fermi-energy'
         RELEN=RELEN-EFERMI
         IF(ABS(EFERMI-EFTEST2)>0.1) PRINT*,'(hope you know what',
     &      ' you are doing!)'
       ENDIF

C--write headers
       PRINT*,'I got for the fermi-energy the value ',EFTEST2

       WRITE(12,*) 'Kriterium 1 (%1): ',(ION(1,I),I=1,ION(1,0)),
     &             ' > ',(CRI(1,I),I=1,3)
       WRITE(12,*) 'Kriterium 2 (%2): ',(ION(2,I),I=1,ION(2,0)),
     &             ' > ',(CRI(2,I),I=1,3)
       IF (EFERMI.LT.-1e9) THEN
         WRITE(12,*) 'calculated and used E(fermi) = ',EFTEST2
       ELSE
         WRITE(12,*) 'given and used E(fermi) = ',EFERMI,
     &               ' (calculated ',EFTEST2,')'
       ENDIF
       WRITE(12,*)
       WRITE(12,*) 'k-Pkt   E     Band  %1  %2  N-8 N-7 N-6 N-5',
     &             ' N-4 N-3 N-2 N-1  N   tot'

C--now calculate surface states and write them out
       DO 2400 NK=1,NKPTS
       DO 2400 NB=1,NBANDS

C I think it is nessecary to calculate PTOT
         PTOT=0.
         DO 2301 I=1,NIONS
2301	   PTOT=PTOT+PARSUM(NK,NB,I)

         DO I=1,2
           PSUM(I)=0
           DO II=1,ION(I,0)
             PSUM(I)=PSUM(I)+PARSUM(NK,NB,ION(I,II))
           ENDDO
           PSUM(I)=PSUM(I)/PTOT
         ENDDO

         IF(PSUM(1).GT.CRI(1,1)) THEN
           WRITE(11,*) 3,NK,RELEN(NK,NB)
           ELSE IF (PSUM(1).GT.CRI(1,2)) THEN
             WRITE(11,*) 2,NK,RELEN(NK,NB)
           ELSE IF (PSUM(1).GT.CRI(1,3)) THEN
             WRITE(11,*) 1,NK,RELEN(NK,NB)
         ENDIF
         IF(PSUM(2).GT.CRI(2,1)) THEN
           WRITE(13,*) 3,NK,RELEN(NK,NB)
           ELSE IF (PSUM(2).GT.CRI(2,2)) THEN
             WRITE(13,*) 2,NK,RELEN(NK,NB)
           ELSE IF (PSUM(2).GT.CRI(2,3)) THEN
             WRITE(13,*) 1,NK,RELEN(NK,NB)
         ENDIF

         IF((PSUM(1).GT.CRI(1,3)).OR.(PSUM(2).GT.CRI(2,3))) THEN
           WRITE(12,'(I3,F10.3,I5,X,3F4.2,X,8F4.2,F5.2)')
     &       NK,RELEN(NK,NB),NB,PSUM(1),PSUM(2),
     &       (PARSUM(NK,NB,I),I=NIONS-8,NIONS),PTOT
         END IF

2400   CONTINUE 	   

2401   CLOSE(10)
       CLOSE(11)
       CLOSE(12)
       CLOSE(13)

      END

C***********************************************************************
C FUNCTION NROFI
C
C gives back the number of integers in the string A
C***********************************************************************
      FUNCTION NROFI(A)
      CHARACTER*256 A
      INTEGER NROFDATA,I,DUM,N

      N=0
      DO
        READ(A,*,ERR=100,END=100) (I,DUM=0,N)
        N=N+1
      ENDDO
 100  NROFI=N
      RETURN
      END

