C-----------------------------------------------------------------------
C File: zeo.for
C Date: 09/12/1997
C-----------------------------------------------------------------------
C Version: 97A
C-----------------------------------------------------------------------
C Tibor F. Nagy, Department of Physics, MSU
C E-mail: nagy_t@pa.msu.edu
C-----------------------------------------------------------------------
C THE PROGRAM:
C - creates "POSITIVE" atomic framework picture
C STEPS:
C - reads atomic radii from "atom.rad"
C - reads the structure from Brookhaven-data-file (extension: ".bkv")
C - calculates the field 
C - writes the field to an AVS ".dat"-file
C-----------------------------------------------------------------------
      PROGRAM     ZEO
C-----------------------------------------------------------------------
      IMPLICIT    NONE
      REAL*4      FIELD(101,101,101)   ! AVS volume field
      REAL*4      XYZ(19999,3)         ! coordinates of the atoms
      REAL*4      RAD(19999)           ! radii of the atoms
      INTEGER     NATOM                ! num of atoms (max: 19999)
      INTEGER     NXGRID               ! number of gridpoints IN X,
      INTEGER     NYGRID               ! ...Y,
      INTEGER     NZGRID               ! and Z directions
      REAL*4      ATOMRAD(99)          ! atomic vdW radii-table
      INTEGER     ATOMNUM(99)          ! atomic number-table
      CHARACTER*2 ATOMNAME(99)         ! atomic name-table
      REAL*4      XMIN,XMAX            ! minimal and maximal...
      REAL*4      YMIN,YMAX            ! positions of the structure...
      REAL*4      ZMIN,ZMAX            ! in the real space
      CHARACTER   FILNAM*30            ! filename
      INTEGER     LENG                 ! length of the filename
C-----------------------------------------------------------------------
C Read the periodic table: atomic radius, number and name
      CALL RDATOMRAD (ATOMRAD,ATOMNUM,ATOMNAME)
C Read the name of the inputfile
      WRITE(6,'(1X,A)') '--- Inputfile: Brookhaven data file ---'
      CALL INQNAM (FILNAM,LENG)
C Read the structure from the Brookhaven data-file
      CALL RDFILE (FILNAM,LENG,ATOMRAD,ATOMNAME,NATOM,XYZ,RAD)
C Put the structure into a cubic grid
      CALL AUTOSC (XYZ,NATOM,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,
     #             NXGRID,NYGRID,NZGRID)
C Calculate the AVS-field
      CALL CLCFLD (NATOM,XYZ,RAD,NXGRID,NYGRID,NZGRID,
     #             XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,FIELD)
C Read the name of the outputfile
      WRITE(6,'(1X,A)') '--- Outputfile: AVS ".dat" file ---'
      CALL INQNAM (FILNAM,LENG)
C Write the field to an AVS ".dat"-file
      CALL WRFILE (FILNAM,LENG,NXGRID,NYGRID,NZGRID,FIELD)
      STOP 
      END
C End of program ZEO
C-----------------------------------------------------------------------
      SUBROUTINE RDATOMRAD (ATOMRAD,ATOMNUM,ATOMNAME)
C-----------------------------------------------------------------------
C Routine to read atomic radii from the file 'atom.rad'
C-----------------------------------------------------------------------
      IMPLICIT    NONE
      REAL*4      ATOMRAD(99)
      INTEGER     ATOMNUM(99)
      CHARACTER*2 ATOMNAME(99)
C Local:
      INTEGER     LIN,I
C-----------------------------------------------------------------------
      LIN = 8
      I = 1
      OPEN(LIN,FILE='atom.rad',STATUS='OLD')
   5  READ(LIN,'(F6.3,1X,I2,1X,A2)',ERR=10,END=10) ATOMRAD(I),
     #ATOMNUM(I),ATOMNAME(I)
C Check the integrity of 'atom.rad' file.
      IF (I.NE.ATOMNUM(I)) THEN
        WRITE(6,'(1X,A)') 'Wrong atomic number in the "atom.rad" file!'
        WRITE(6,'(1X,A,I4)') 'Line number: ',I
        STOP
      ENDIF
      I = I + 1
      GOTO 5  !go back to read next one
  10  IF (I.NE.100) THEN
        WRITE (6,'(1X,A)') 'The "atom.rad" file is not complete.'
        WRITE (6,'(1X,A,I4)') 'Line number: ',I
        STOP
      ENDIF
      CLOSE(LIN)
      RETURN
      END
C End of subroutine RDATOMRAD
C-----------------------------------------------------------------------
      SUBROUTINE INQNAM (FILNAM,LENG)
C-----------------------------------------------------------------------
C Read a string and compute the length of the string
C-----------------------------------------------------------------------
      IMPLICIT   NONE
      CHARACTER  FILNAM*30
      INTEGER    LENG
C-----------------------------------------------------------------------
   1  WRITE(6,'(1X,A,$)') 'Enter filename. (no extension): '
      READ (5,'(A30)') FILNAM
      LENG = 30
      DO WHILE (FILNAM(LENG:LENG).EQ.' ')
        LENG = LENG - 1
      ENDDO
      IF (LENG.LT.1) GOTO 1
      RETURN
      END
C End of subroutine INQNAM
C-----------------------------------------------------------------------
      SUBROUTINE RDFILE (FILNAM,LENG,ATOMRAD,ATOMNAME,NATOM,XYZ,RAD)
C-----------------------------------------------------------------------
C Routine to read Brookhaven file 
C (Clean and simple version.)
C-----------------------------------------------------------------------
      IMPLICIT     NONE
      CHARACTER*30 FILNAM
      INTEGER      LENG
      REAL*4       ATOMRAD(99)
      CHARACTER*2  ATOMNAME(99)
      INTEGER      NATOM
      REAL*4       XYZ(19999,3),RAD(19999)
C Local:
      INTEGER      LIN             ! logical I/O channel
      REAL         X,Y,Z,RADIUS    ! coordinates and radius 
C Local variables for Brookhaven Database
      CHARACTER    SPEC*6          ! start string
      CHARACTER    INPUT*70        ! input string 
      CHARACTER    DUM1*6,DUM2*15  ! dummy variables for unnecessary fields
      CHARACTER*1  C1,C2,C3        ! three characters for atomname
      LOGICAL      F1,F2,F3        ! three logical variables
      INTEGER      I1,I2,I3        ! three Ascii codes
      CHARACTER*2  CNAME           ! current atom name
      INTEGER      I               ! loop variable
C-----------------------------------------------------------------------
C Logical input for ".bkv"-file
      LIN = 8
C Reset atom-counter
      NATOM = 0
C Message to the screen
      WRITE(6,'(1X,A)') 'Reading atoms from BKV file...'
C Open Brookhaven file
      OPEN(UNIT=LIN,FILE=FILNAM(1:LENG)//'.bkv',STATUS='OLD')
  5   READ(LIN,'(A6,A70)',ERR=10) SPEC,INPUT
      IF (((SPEC(1:6).EQ.'HEADER').OR.(SPEC(1:6).EQ.'REMARK'))) GOTO 5
      IF (((SPEC(1:6).EQ.'CONECT').OR.(SPEC(1:6).EQ.'END   '))) GOTO 10
      IF (((SPEC(1:6).EQ.'HETATM').OR.(SPEC(1:6).EQ.'ATOM  '))) THEN
        READ(INPUT,99) DUM1,C1,C2,C3,DUM2,X,Y,Z
 99     FORMAT(A6,3A1,A15,3F8.3)
C Testing the three characters
        CALL LETTERP(C1,F1,I1)
        CALL LETTERP(C2,F2,I2)
        CALL LETTERP(C3,F3,I3)
C Building the elementname from the three characters
        IF (F1) THEN
          IF (F2) THEN
            CNAME = C1//C2
          ELSE
            CNAME = C1//' '
          ENDIF
        ELSE
          IF (F2) THEN
            IF (F3) THEN
              CNAME = C2//C3
            ELSE
              CNAME = C2//' '
            ENDIF
          ELSE
            WRITE(6,'(1X,A)') 'Illegal atomname-field in BKV-file.'
            WRITE(6,'(1X,A,I6)')  'Natom: ',(NATOM+1)
            WRITE(6,'(1X,A,3A1)') 'Field: ',C1,C2,C3
            STOP
          ENDIF
        ENDIF
C Elementname identification.
        DO I=1,99
          IF (ATOMNAME(I).EQ.CNAME) THEN
            RADIUS=ATOMRAD(I)
            GOTO 25
          ENDIF
        ENDDO
        WRITE(6,'(1X,A)')    'Error: Atom type is unknown!'
        WRITE(6,'(1X,2A)')   'Name: ',CNAME
        WRITE(6,'(1X,A,I6)') 'Atom: ',(NATOM+1)
        STOP
C The name was found to be valid. Let's store the atom.
 25     NATOM = NATOM + 1
        IF(NATOM.GT.19999) THEN
          WRITE(6,'(1X,A)') 'Error: too many atoms'
          STOP
        ENDIF
C Let's store the coordinates and radius
        XYZ(NATOM,1) = X
        XYZ(NATOM,2) = Y
        XYZ(NATOM,3) = Z
        RAD(NATOM)   = RADIUS
C Test module.
C        write(6,*) 'natom: ',natom
C        write(6,*) 'xyzr : ',x,y,z,radius
C        write(6,*) 'name : ',cname
C Test module ends.
        GOTO 5     ! Go back to read in next atom
      ENDIF
 10   CLOSE(LIN)   ! end of file reached, finish.
      WRITE(6,'(1X,A,I6)') 'Number of atoms in structure: ',NATOM
      RETURN
      END
C End of subroutine RDFILE
C-----------------------------------------------------------------------
      SUBROUTINE LETTERP(LETTER,FLAG,CODE)
C-----------------------------------------------------------------------
C Subroutine to test a character if it is a capital or small letter.
C If it is a small letter, convert it to capital.
C If it is not a letter at all, the FLAG will be .FALSE.
C-----------------------------------------------------------------------
      IMPLICIT  NONE
      CHARACTER LETTER*1 ! Input character to be tested
      LOGICAL   FLAG     ! TRUE if LETTER is a capital or small letter
      INTEGER   CODE     ! Ascii code of the LETTER
C-----------------------------------------------------------------------
      CODE = ICHAR(LETTER)
      IF ((CODE.GE.97).AND.(CODE.LE.122)) THEN 
        CODE = CODE - 32
        LETTER = CHAR(CODE)
      ENDIF
      FLAG = ((CODE.GE.65).AND.(CODE.LE.90))
      RETURN
      END
C End of subroutine LETTERP
C-----------------------------------------------------------------------
      SUBROUTINE AUTOSC (XYZ,NATOM,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,
     #                   NXGRID,NYGRID,NZGRID)
C-----------------------------------------------------------------------
C Put the structure or molecule in a cubic shape grid
C-----------------------------------------------------------------------
      IMPLICIT  NONE
      REAL*4    XYZ(19999,3)
      INTEGER   NATOM
      INTEGER   NXGRID,NYGRID,NZGRID
      REAL*4    XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
C Local:
      REAL*4    MIDX,MIDY,MIDZ
      REAL*4    DELX,DELY,DELZ,DEL,EDGE
      REAL*4    EXFACT
      INTEGER   I
C-----------------------------------------------------------------------
C Read AVS resolution
 10   WRITE(6,'(1X,A22,A30,$)') 'Enter AVS resolution. ',
     #'(min:1 / max:101 / optimal:41) '
      READ(5,*) NXGRID
      IF ((NXGRID.GT.101).OR.(NXGRID.LT.1)) GOTO 10
      NYGRID = NXGRID
      NZGRID = NXGRID
C Read expanding factor
      WRITE(6,'(1X,A24,A37,$)') 'Enter Expanding Factor. ',
     #'(1:normal / LT1:magnif / GT1:reduce) '
      READ(5,*) EXFACT
      XMIN = XYZ(1,1)
      XMAX = XYZ(1,1)
      YMIN = XYZ(1,2)
      YMAX = XYZ(1,2)
      ZMIN = XYZ(1,3)
      ZMAX = XYZ(1,3)
      DO I=2,NATOM
        IF (XYZ(I,1).LT.XMIN) XMIN=XYZ(I,1)
        IF (XYZ(I,1).GT.XMAX) XMAX=XYZ(I,1)
        IF (XYZ(I,2).LT.YMIN) YMIN=XYZ(I,2)
        IF (XYZ(I,2).GT.YMAX) YMAX=XYZ(I,2)
        IF (XYZ(I,3).LT.ZMIN) ZMIN=XYZ(I,3)
        IF (XYZ(I,3).GT.ZMAX) ZMAX=XYZ(I,3)
      ENDDO
C Middle points
      MIDX = (XMAX+XMIN) / 2.0
      MIDY = (YMAX+YMIN) / 2.0
      MIDZ = (ZMAX+ZMIN) / 2.0
C Sizes in different directions
      DELX = XMAX - XMIN
      DELY = YMAX - YMIN
      DELZ = ZMAX - ZMIN
C Biggest dimension
      DEL  = MAX(DELX,DELY,DELZ) / 2.0
      XMIN = MIDX - DEL*EXFACT
      XMAX = MIDX + DEL*EXFACT
      YMIN = MIDY - DEL*EXFACT
      YMAX = MIDY + DEL*EXFACT
      ZMIN = MIDZ - DEL*EXFACT
      ZMAX = MIDZ + DEL*EXFACT
      EDGE = 2.0*DEL*EXFACT
      WRITE(6,'(1X,A,F12.3,A)') 'Cube-Edge-Length : ',EDGE,' A'
      WRITE(6,'(1X,A,F12.3,A)') 'Total Cube Volume: ',EDGE**3.,' A^3'
      RETURN
      END
C End of subroutine AUTOSC
C-----------------------------------------------------------------------
      SUBROUTINE CLCFLD (NATOM,XYZ,RAD,NXGRID,NYGRID,NZGRID,
     #                   XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,FIELD)
C-----------------------------------------------------------------------
C Routine to calculate the scalar-field
C-----------------------------------------------------------------------
      IMPLICIT    NONE
      REAL*4      XYZ(19999,3),RAD(19999)
      INTEGER     NATOM
      INTEGER     NXGRID,NYGRID,NZGRID
      REAL*4      XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
      REAL*4      FIELD(101,101,101)
C Local:
      REAL*4      RAD2(19999)
      INTEGER     I,J,K,L,COUNT
      REAL*4      XSTEP,YSTEP,ZSTEP
      REAL*4      XP,YP,ZP
      REAL*4      XW,YW,ZW,RW
      REAL*4      RATIO
C-----------------------------------------------------------------------
C Reset the we-are-inside-of-the-atom counter 
      COUNT = 0
C Square the atomic radii
      DO I=1,NATOM
        RAD2(I) = RAD(I)*RAD(I)
      ENDDO
C Step sizes
      XSTEP = (XMAX-XMIN)/(NXGRID-1)
      YSTEP = (YMAX-YMIN)/(NYGRID-1)
      ZSTEP = (ZMAX-ZMIN)/(NZGRID-1)
      WRITE(6,'(1X,A,I2)') 'Counting to ',NXGRID
      WRITE(6,'(1X,I2)') 0
C Three loops over the three directions
      DO I=1,NXGRID
        IF (MOD(I,5).EQ.0) WRITE(6,'(1X,I2)') I
        XP = XMIN + (I-1)*XSTEP
        DO J=1,NYGRID
          YP = YMIN + (J-1)*YSTEP
          DO K=1,NZGRID
            ZP = ZMIN + (K-1)*ZSTEP
            FIELD(I,J,K) = 0.0
C Loop over the atoms
            DO L=1,NATOM
              XW = XP - XYZ(L,1)
              YW = YP - XYZ(L,2)
              ZW = ZP - XYZ(L,3)
              RW = XW*XW + YW*YW + ZW*ZW
              IF (RW.LT.RAD2(L)) THEN
                FIELD(I,J,K) = 1.0
                COUNT = COUNT + 1
                GOTO 10
              ENDIF
            ENDDO
 10       ENDDO
        ENDDO
      ENDDO
      RATIO = FLOAT(COUNT) / FLOAT(NXGRID*NYGRID*NZGRID)
      WRITE(6,'(1X,A)') 'Ratios:'
      WRITE(6,'(1X,A,F6.3)') 'Atom/Total : ',RATIO
      WRITE(6,'(1X,A,F6.3)') 'Empty/Total: ',(1.0-RATIO)
      RETURN
      END
C End of subroutine CLCFLD
C-----------------------------------------------------------------------
      SUBROUTINE WRFILE (FILNAM,LENG,NXGRID,NYGRID,NZGRID,FIELD)
C-----------------------------------------------------------------------
C Routine to write field to an AVS/Explorer ".dat"-file
C-----------------------------------------------------------------------
      IMPLICIT   NONE
      REAL*4     FIELD(101,101,101)
      INTEGER    NXGRID,NYGRID,NZGRID
      CHARACTER  FILNAM*30
      INTEGER    LENG
C Function to put one character to a file (Unix only!)
      INTEGER    FPUTC
C Used functions
      REAL*4     FLDMIN,FLDMAX
C Local:
      INTEGER    LOUT,I,J,K
      INTEGER    IER,IW
      REAL*4     FMIN,FMAX
      REAL*4     SCALE,FW
C Dummy array for AVS file
      REAL*4     ZERO(13)
      DATA       ZERO /13*0.0/
C-----------------------------------------------------------------------
C Delete the outest surface for visualization
      DO I=1,NXGRID
        DO J=1,NYGRID
          FIELD(I,J,1)      = 0.0
          FIELD(I,J,NZGRID) = 0.0
        ENDDO
      ENDDO
      DO I=1,NXGRID
        DO J=1,NZGRID
          FIELD(I,1,J)      = 0.0
          FIELD(I,NYGRID,J) = 0.0
        ENDDO
      ENDDO
      DO I=1,NYGRID
        DO J=1,NZGRID
          FIELD(1,I,J)      = 0.0
          FIELD(NXGRID,I,J) = 0.0
        ENDDO
      ENDDO
C Search the min and max value of the field for scaling
      FMIN = FLDMIN(FIELD,NXGRID,NYGRID,NZGRID)
      FMAX = FLDMAX(FIELD,NXGRID,NYGRID,NZGRID)
      WRITE(6,'(1X,A)') 'Scaling:'
      WRITE(6,'(1X,A,F12.3,A)') 'Maximal field-value :',FMAX,
     #' : 255'
      WRITE(6,'(1X,A,F12.3,A)') 'Minimal field-value :',FMIN,
     #' :   0'
C Scaling factor
      SCALE = 255.0 / (FMAX-FMIN)
C Logical output unit
      LOUT = 8
      OPEN(LOUT,FILE=FILNAM(1:LENG)//'.dat',FORM='UNFORMATTED',
     #     STATUS='UNKNOWN')
C 1st: write field dimensions to the file
      IER = FPUTC(LOUT,CHAR(NXGRID))
      IER = FPUTC(LOUT,CHAR(NYGRID))
      IER = FPUTC(LOUT,CHAR(NZGRID))
C 2nd: write the field
      DO I=1,NZGRID
         DO J=1,NYGRID
            DO K=1,NXGRID
               FW  = FIELD(K,J,I)
               IW  = SCALE*(FW-FMIN)
               IER = FPUTC(LOUT,CHAR(IW))
            ENDDO
         ENDDO
      ENDDO
C 3rd: write some more parameters for AVS
      WRITE(LOUT) 0,ZERO,0,FMIN,255,FMAX
      CLOSE(LOUT)
      RETURN
      END
C End of subroutine WRFILE
C-----------------------------------------------------------------------
      REAL*4 FUNCTION FLDMIN(ARRAY,NX,NY,NZ)
C-----------------------------------------------------------------------
C Function to calculate the minimum value in a 3d-array
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*4   ARRAY(101,101,101)
      INTEGER  NX,NY,NZ
C Local:
      INTEGER  I,J,K
C-----------------------------------------------------------------------
      FLDMIN = ARRAY(1,1,1)
      DO I=1,NZ
        DO J=1,NY
          DO K=1,NX
            IF (ARRAY(K,J,I).LT.FLDMIN) FLDMIN=ARRAY(K,J,I)
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END
C End of function FLDMIN
C-----------------------------------------------------------------------
      REAL*4 FUNCTION FLDMAX(ARRAY,NX,NY,NZ)
C-----------------------------------------------------------------------
C Function to calculate the maximum value in a 3d-array
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*4   ARRAY(101,101,101)
      INTEGER  NX,NY,NZ
C Local:
      INTEGER  I,J,K
C-----------------------------------------------------------------------
      FLDMAX = ARRAY(1,1,1)
      DO I=1,NZ
        DO J=1,NY
          DO K=1,NX
            IF (ARRAY(K,J,I).GT.FLDMAX) FLDMAX=ARRAY(K,J,I)
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END
C End of function FLDMAX
C-----------------------------------------------------------------------
C Je me souviens
C Des jours anciens
C Et je pleure.
C-----------------------------------------------------------------------
C23456789012345678901234567890123456789012345678901234567890123456789012
C-----------------------------------------------------------------------