      Program chgprp                                                    CHGP001
c --------------------------------------------------------------------- CHGP002
c         prepares charge.inp input file from SEXIE output files        CHGP003
c                    CORMAT.DAT, CELL.DAT, INFEFF.XYZ                   CHGP004
c      infeff.xyz is actually written in the subroutine priver.for      CHGP005
c               see priver.for for revisions to infeff.xyz              CHGP006
c --------------------------------------------------------------------- CHGP007
c     VARIABLES:                                                        CHGP008
c       nat   =  number of atoms in shell                               CHGP009
c       x,y,z =  cartesion coordinates                                  CHGP010
c       charge=  charge of ligand atom                                  CHGP011
c       foc   =  occupation factor                                      CHGP012
c       xnn   =  distance from central atom                             CHGP013
c       el    =  element name                                           CHGP014
c       iele  =  atomic number                                          CHGP015
c       rmax  =  largest permitted radius                               CHGP016
c       iweigt=  weighting 1/r**iweigt                                  CHGP017
c     Compatible with SEXIE 5.0 - B.Rupp (1995)                         CHGP018
c --------------------------------------------------------------------- CHGP019
      parameter (iat=20, ims=100, maxats=48)                            CHGP020
c --- maximum number of atoms in shell is maxats = 48   (SEXIE limit)   CHGP021
c     maximum number of shells to be read in ims = 100  (SEXIE limit)   CHGP022
c     maximum number of inequivalent atoms iat   = 20   (SEXIE limit)   CHGP023
c     NOTE : SEXIE limits need to be changed in all subroutines !       CHGP024
c     (limits should be higher for low symmetry problems)               CHGP025
                                                                        CHGP026
      character inname*10,outname*10,el*2,revdat*8,icomp*4,ipotid*4     CHGP027
      character doc*60,dc1*60,dc2*60,title*60,titlesg*60,ipotat*2,ide*4 CHGP028
      logical exists                                                    CHGP029
      dimension x(ims+1,maxats), y(ims+1,maxats), z(ims+1,maxats)       CHGP030
      dimension ipot(ims+1),ipotat(iat+1),ipotid(iat+1),ipotel(iat+1)   CHGP031
      dimension nat(ims+1),xnn(ims+1),ide(ims+1),ielec(ims+1),el(ims+1) CHGP032
      dimension charge(ims+1), t(3,3),cel(6),foc(ims+1)                 CHGP033
      inname='infeff.xyz'                                               CHGP034
      outname='charge.inp'                                              CHGP035
                                                                        CHGP036
c --- revision comments and date here ----                              CHGP037
c     version  5.0  by B.Rupp                                           CHGP038
      revdat='08/21/95'                                                 CHGP039
                                                                        CHGP040
c --- look for SEXIE configuration file and read the computer type --   CHGP041
      inquire (file='SEXIE.CFG',exist=exists)                           CHGP042
      if (exists) then                                                  CHGP043
         open (4,file='SEXIE.CFG',status='old')                         CHGP044
         do 0001 i=1,7                                                  CHGP045
 0001       read (4,'(8x,a4)') icomp                                    CHGP046
         read(4,'(8x,a4)')   icomp                                      CHGP047
         call upstrg (icomp,4)                                          CHGP048
         close(4,status='keep')                                         CHGP049
      else                                                              CHGP050
         write (*,'(/a)') '  Cannot find configuration file SEXIE.CFG'  CHGP051
         write (*,*) ' Resume with 7-bit TTY terminal character set'    CHGP052
         icomp  = '7BIT'                                                CHGP053
      end if                                                            CHGP054
                                                                        CHGP055
c --- Initalize ----                                                    CHGP056
      do 0002 i=1,ims+1                                                 CHGP057
          ipot(i)=0                                                     CHGP058
          ielec(i)=0                                                    CHGP059
          xnn(i)=0.                                                     CHGP060
          nat(i)=0                                                      CHGP061
          foc(i)=0.                                                     CHGP062
          if(i.le.(iat+1)) then                                         CHGP063
            ipotat(i)='  '                                              CHGP064
            ipotid(i)='    '                                            CHGP065
          end if                                                        CHGP066
          do 0003 j=1,maxats                                            CHGP067
            x(i,j)=0.                                                   CHGP068
            y(i,j)=0.                                                   CHGP069
            z(i,j)=0.                                                   CHGP070
 0003     continue                                                      CHGP071
 0002 continue                                                          CHGP072
                                                                        CHGP073
c --- header ------                                                     CHGP074
                                                                        CHGP075
      call chghed (revdat,icomp)                                        CHGP076
                                                                        CHGP077
c --- read file infeff.xyz ---------------                              CHGP078
      open (5,file=inname,status='old')                                 CHGP079
      read (5,'(a60)') dc1                                              CHGP080
      read (5,'(a60)') title                                            CHGP081
      read (5,'(a60)') titlesg                                          CHGP082
c --- loop over all shells up to ims +1 ---                             CHGP083
      do 0004 i=1,ims+1                                                 CHGP084
         read (5,'(a60)',end=9001,err=9009) doc                         CHGP085
         read (5,'(a60)',end=9001,err=9009) dc2                         CHGP086
c ---    read how many atoms in this shell ----                         CHGP087
         read (5,*) nat(i)                                              CHGP088
c ---    loop over all atoms in this shell ----                         CHGP089
         do 0005 j=1,nat(i)                                             CHGP090
            read (5,*)x(i,j)                                            CHGP091
            read (5,*)y(i,j)                                            CHGP092
            read (5,*)z(i,j)                                            CHGP093
 0005       continue                                                    CHGP094
         read (5,*) xnn(i)                                              CHGP095
         read (5,'(1x,a2,1x,a4,2f8.4)') el(i),ide(i),charge(i),foc(i)   CHGP096
         read (5,*) ielec(i)                                            CHGP097
                                                                        CHGP098
c ---    sort the atom kinds for potential ----                         CHGP099
c        if element and crystallographic position equivalent, then      CHGP100
c        assign the same potential flag                                 CHGP101
         ihits=0                                                        CHGP102
         if (i.le.2) then                                               CHGP103
c           put first and second shell in potential list                CHGP104
            ipot(i)=i                                                   CHGP105
            ipotat(i)=el(i)                                             CHGP106
            ipotid(i)=ide(i)                                            CHGP107
            ipotel(i)=ielec(i)                                          CHGP108
            ipots=2                                                     CHGP109
         else                                                           CHGP110
c           check if we have new atom already in potential list         CHGP111
            do 0006 ii=2,ipots                                          CHGP112
               if((el(i).eq.ipotat(ii)).and.(ide(i).eq.ipotid(ii)))then CHGP113
                 ihits=1                                                CHGP114
                 ipot(i)=ii                                             CHGP115
               end if                                                   CHGP116
 0006      continue                                                     CHGP117
           if (ihits.eq.0) then                                         CHGP118
              ipots=ipots+1                                             CHGP119
              ipot(i)=ipots                                             CHGP120
              ipotat(ipots)=el(i)                                       CHGP121
              ipotid(ipots)=ide(i)                                      CHGP122
              ipotel(ipots)=ielec(i)                                    CHGP123
           end if                                                       CHGP124
         end if                                                         CHGP125
 0004 continue                                                          CHGP126
 9001 close (5)                                                         CHGP127
c --- decrement loop counter ---                                        CHGP128
      imax=i-1                                                          CHGP129
                                                                        CHGP130
c --- get cell constants --                                             CHGP131
      open(4,file='cell.dat',status='old')                              CHGP132
      read(4,*) (cel(i),i=1,6)                                          CHGP133
      close(4)                                                          CHGP134
c --- get coordinate transformation matrix --                           CHGP135
      open(3,file='cormat.dat',status='old')                            CHGP136
      do 0007 j=1,3                                                     CHGP137
         read(3,'(3f12.6)') (t(j,i),i=1,3)                              CHGP138
 0007 continue                                                          CHGP139
      close(3)                                                          CHGP140
                                                                        CHGP141
c --- screen output ----  
      write (*,*)                                                       CHGP142
      write (*,'(a,i3)')'  Number of shells read : ', imax-1            CHGP143
      write (*,'(a,f7.4)') '  Radius of last shell read : ',xnn(imax-1) CHGP144
                                                                        CHGP145
c --- write output file                                                 CHGP146
      open (6,file=outname, status='unknown')                           CHGP147
      write (6,*)'****************************************************' CHGP148
      write (6,*)'*     CHARGE input file generated by SEXIE         *' CHGP149
      write (6,*)'****************************************************' CHGP150
      write (6,*)'TITLE  ',title                                        CHGP151
      write (6,*)'TITLE  ',titlesg                                      CHGP152
      write (6,'(1x,a7,3f8.4,3f7.2)')'CELL   ',(cel(i), i=1,6)          CHGP153
      do 0010 j=1,3                                                     CHGP154
         write (6,'(1x,a7,3f12.6)')'CORMAT ',  (t(j,i),i=1,3)           CHGP155
 0010 continue                                                          CHGP156
      write (6,'(1x,a7,10i3)')'SHELLS ',0,0,0,0,0,0,0,0,0,0             CHGP157
      write (6,'(1x,a7,f8.4,a)')'RMAX   ',0.0,'   Max. radial distance' CHGP158
      write (6,'(1x,a7,i8,a)')'WEIGHT ',0,'   Weight (1/r**n)'          CHGP159
      write (6,'(1x,a7,i8,a)')'USECHG ',0,'   Use charge (1) do not (0)'CHGP160
c --- atom list ---                                                     CHGP161
      write (6,*)'CENTRE'                                               CHGP162
      i=1                                                               CHGP163
      j=1                                                               CHGP164
      write (6,'(i4,1x,a4,a4,6f8.4)') i-1,el(i),                        CHGP165
     &      ide(i),x(i,j),y(i,j),z(i,j),xnn(i),charge(i),foc(i)         CHGP166
      write (6,'(a,/a)')' ATOMS ',                                      CHGP167
     &' *Shl atm site    x       y       z    dist   charge    foc'     CHGP168
c --- loop over all shells                                              CHGP169
      do 0008 i=2,imax                                                  CHGP170
c --- and print all atoms of the shell                                  CHGP171
        do 0009 j=1,nat(i)                                              CHGP172
           write (6,'(i4,1x,a4,a4,6f8.4)') i-1,el(i),                   CHGP173
     &      ide(i),x(i,j),y(i,j),z(i,j),xnn(i),charge(i),foc(i)         CHGP174
 0009   continue                                                        CHGP175
 0008 continue                                                          CHGP176
      write (6,*) 'END'                                                 CHGP177
      close (6)                                                         CHGP178
      stop ' Program terminated ok'                                     CHGP179
 9009 write (*,'(/a,a,a/)') ' Error reading ', inname,' - please check' CHGP180
      end                                                               CHGP181
