                          PROGRAM SEXIE                                 SEXI001
c ----------------------------------------------------------------------SEXI002
c     This program module is part of the PhysiSoft X-ray Library        SEXI003
c     and may not be sold or reproduced for commercial purposes.        SEXI004
c     Although it has been tested carefully, there is no warranty       SEXI005
c     regarding accuracy or freedom from side-effects.                  SEXI006
c                   (C) 1986-1994 Bernhard Rupp                         SEXI007
c ----------------------------------------------------------------------SEXI008
c                      Generic CPC FORTRAN version                      SEXI009
c ----------------------------------------------------------------------SEXI010
c  A. Input file : SEXIE.INP (unit=5)                                   SEXI011
c                                                                       SEXI012
c    LINE 1: A(68)                                                      SEXI013
c     title  : title line , 68 characters                               SEXI014
c    LINE 2: F,F,I,F                                                    SEXI015
c     dmax   : maximum distance up to which coord.cluster is calculated SEXI016
c     condis : range within which atoms are considered to               SEXI017
c              be in same shell                                         SEXI018
c     isearc : search flag for polyhedra determination                  SEXI019
c     conang : range within angles are considered to be special angles  SEXI020
c    LINE 3: 6(F)                                                       SEXI021
c     cel(1-6) : a,b,c,al,be,ga : cell constants                        SEXI022
c    LINE 4: A(10,I4)                                                   SEXI023
c     spcgrp : space group symbol and number (see INSHELL for details)  SEXI024
c    LINES 5 to EOF : fixed (A4,A4,3f8.5,I3)                            SEXI025
c     elemt  : atom name (iat)                                          SEXI026
c     ide    : atom label (iat)                                         SEXI027
c     x,y,z  : atom positions (iat)                                     SEXI028
c     icalc  : flag for calculation of shells for corresp. atom (iat)   SEXI029
c                                                                       SEXI030
c  B. Variables used in program :                                       SEXI031
c     iatom  : number of inequivalent atom positions  (max iat)         SEXI032
c     ns     : number of point transformations (max iqv)                SEXI033
c     ts     : translation vector of equipoint transformation           SEXI034
c     rs     : rotation martix of equipoint transformation              SEXI035
c     imult  : multiplicity of position (iat)                           SEXI036
c     igen   : multiplicity of general position                         SEXI037
c     iortho : flag for orthogonal systems 1                            SEXI038
c     iwarn  : contains warning flag if equipoint atoms coincide        SEXI039
c     st     : indexed array stores all generated atom positions        SEXI040
c              and their distance to origin                             SEXI041
c     ipt    : pointer to st                                            SEXI042
c     icntra : pointer to st                                            SEXI043
c     did    : position, distance, atom code, cart.coord  (istor,8)     SEXI044
c     ishell : coord.number of shell (ish)                              SEXI045
c     nshell : number of coord. shells within d(max)                    SEXI046
c     cosi,sinu,cosq : cosine, sine and cosine**2 of cell angles        SEXI047
c     t      : contravariant transformation matrix into orthonormal     SEXI048
c     isub   : index array for subshell components                      SEXI049
c ----------------------------------------------------------------------SEXI050
c     This program minimizes memory requirements so that it can         SEXI051
c     run on simple PC's as well. There are several pointers to         SEXI052
c     arrays (and pointers to pointers), which are not explicitly       SEXI053
c     described in the documentation. Please take care whenever         SEXI054
c     changing assignment of variables which might be indexing          SEXI055
c     variables (arrays) connected to primary data arrays.              SEXI056
c     For easier modifications, no common blocks are used and all       SEXI057
c     global variables are explicitly referenced by name or value       SEXI058
c ----------------------------------------------------------------------SEXI059
c     Any improvements and critizism (in particular for a fast          SEXI060
c     mathematical algorithm to find the coordination geometry)         SEXI061
c     should be sent to :                                               SEXI062
c     Bernhard Rupp                                                     SEXI063
c     BR@OEDIPUS.LLNL.GOV                                               SEXI064
c     LLNL - BBRP                                                       SEXI065
c ----------------------------------------------------------------------SEXI066
c     Parameters used in main program and most subroutines:             SEXI067
c     iat    : max. number of nonequivalent sites or atoms in cell (20) SEXI068
c     istor  : max. number of atoms in cluster (2500)                   SEXI069
c     ims    : max. number of shells (100)                              SEXI070
c     ipr    : print file unit 6 - do not change                        SEXI071
c ----------------------------------------------------------------------SEXI072
c     files opened : 3 = infeff.xyz, shell info for J.Rehr's FEFF       SEXI073
c                    4 = configuration file SEXIE.CFG                   SEXI074
c                    5 = sexie.inp                                      SEXI075
c                    6 = screen                                         SEXI076
c                    7 = sexie.out, standard output file                SEXI077
c                    8 = xyz.dat and atom.dat, contain data for SCHAKAL SEXI078
c                    9 = cell.dat, contains unit cell data via CELOUT   SEXI079
c ----------------------------------------------------------------------SEXI080
                                                                        SEXI081
      parameter (iat=20,istor=2500,ims=100,ipr=6)                       SEXI082
      character elemt*2,ide*4,title*60,revdat*8,icomp*4,avern*3         SEXI083
      logical exists                                                    SEXI084
      dimension elemt(iat),ide(iat),x(iat),y(iat),z(iat),igr(36)        SEXI085
      dimension ts(3,24),fs(3,3,24),imult(iat)                          SEXI086
      dimension iwarn(iat),st(iat,192,4),ipt(iat),icalc(iat)            SEXI087
      dimension icntra(iat),did(istor,8),ish(ims),dist(ims)             SEXI088
      dimension cel(6),cosi(3),sinu(3),t(3,3),cosq(3)                   SEXI089
      dimension isub(ims,iat)                                           SEXI090
                                                                        SEXI091
c --- revision data -- PUT ANY REVISION MARKS IN HERE PLEASE ---------  SEXI092
      revdat='02/10/94'                                                 SEXI093
      avern='3.0'                                                       SEXI094
c     last revision includes Subroutine CELOUT, PRIVER :                SEXI095
c     numerical error in transformation from fractional to cartesian    SEXI096
c     coordinates - fixed                                               SEXI097
c     basis tramsformation made standard with IUCr                      SEXI098
c     I41/amd and P3221 decoded wrong due to rounding precision-fixed   SEXI099
c     all system files upper case now                                   SEXI100
c     PRIVER.FOR reads ispa now -  A.T.M.                               SEXI101
                                                                        SEXI102
c --- look for SEXIE configuration file and read the computer type --   SEXI103
      inquire (file='SEXIE.CFG',exist=exists)                           SEXI104
      if (exists) then                                                  SEXI105
         open (4,file='SEXIE.CFG',status='old')                         SEXI106
         do 8808 i=1,7                                                  SEXI107
 8808       read (4,'(8x,a4)') icomp                                    SEXI108
         read(4,'(8x,a4)')   icomp                                      SEXI109
         call upstrg (icomp,4)                                          SEXI110
         close(4,status='keep')                                         SEXI111
      else                                                              SEXI112
         write (*,'(/a)') '  Cannot find configuration file SEXIE.CFG'  SEXI113
         write (*,*) ' Resume with 7-bit TTY terminal character set'    SEXI114
         icomp  = '7BIT'                                                SEXI115
      end if                                                            SEXI116
                                                                        SEXI117
c --- defauft unit screen = 6, output file = 7 -                        SEXI118
      open (unit=ipr+1,file='sexie.out',status='unknown')               SEXI119
                                                                        SEXI120
c --- reset transformation matrix, space group symbol and div. flags    SEXI121
      call reset (ts,fs,igr,isyst,ns,isymce,iwarn,icalc)                SEXI122
                                                                        SEXI123
c --- read input file --                                                SEXI124
      call getinp (title,elemt,ide,x,y,z,cel,iatom,igr,                 SEXI125
     &             isearc,conang,condis,dmax,ispa,icalc)                SEXI126
                                                                        SEXI127
c --- reset if unrealistic parameters for shell search -                SEXI128
      if (condis.lt.0.00001) condis=0.00001                             SEXI129
      if (conang.lt.0.05) conang=0.05                                   SEXI130
                                                                        SEXI131
c --- retrieve system from cell constants -                             SEXI132
      call system (cel,isystm,iortho)                                   SEXI133
                                                                        SEXI134
c --- calculate trigonometric terms for storage -                       SEXI135
      call trigo(cel,cosi,sinu,cosq)                                    SEXI136
                                                                        SEXI137
c --- calculate cell volume  -                                          SEXI138
      call volume (cel,cosq,cosi,vez,iortho)                            SEXI139
                                                                        SEXI140
c --- calculate the transformation matrix T for cell to orthonormal -   SEXI141
      call trigo2 (cosi,sinu,t,cosq)                                    SEXI142
                                                                        SEXI143
c --- decode space group for equipoint transformations                  SEXI144
      call equipt (isyst,isymce,isymbr,ns,ts,fs,igr)                    SEXI145
                                                                        SEXI146
c --- get number of bravais lattice -                                   SEXI147
      ibravl=ibrav(isymbr)                                              SEXI148
                                                                        SEXI149
c --- print first part of output while checking consistency of system   SEXI150
      call oneout (title,dmax,conang,condis,isearc,iatom,vez,cel,       SEXI151
     &    isystm,isyst,igr,isymce,ibravl,ispa,t,cosi,revdat,icomp,avern)SEXI152
                                                                        SEXI153
c --- print equipoint matrix and decoded equipoints -                   SEXI154
      call eqiout (ts,fs,ns,isymce,ibravl,igen,icomp)                   SEXI155
                                                                        SEXI156
c --- end program if no atom card found --                              SEXI157
      if (iatom.eq.0) goto 9001                                         SEXI158
                                                                        SEXI159
c --- multiplicities of atom positions and print atoms in asym.cell --  SEXI160
      call multip (iatom,ibravl,x,y,z,fs,ts,imult,isymce,ns,igen,       SEXI161
     &             iwarn,cel,st,ipt)                                    SEXI162
                                                                        SEXI163
c --- find the central atom for expansion closest to (000) --           SEXI164
      call centra(iatom,ipt,st,cel,cosi,icntra,dismax)                  SEXI165
                                                                        SEXI166
c --- expand the matrix T with the cell dimensions --                   SEXI167
      call trigo3(cel,t)                                                SEXI168
                                                                        SEXI169
c --- print atom positions and unit cell information -                  SEXI170
      call celout (iatom,elemt,ide,x,y,z,imult,iwarn,st,ipt,dismax,     SEXI171
     &             vez,t,cel,icomp,title)                               SEXI172
                                                                        SEXI173
c --- end program here, if dmax=lt 3.0 --                               SEXI174
      if (dmax.lt.3.0) goto 9002                                        SEXI175
                                                                        SEXI176
c --- end program here if no shells to calculate --                     SEXI177
      isum=0                                                            SEXI178
      do 0003 i=1,iat                                                   SEXI179
 0003    isum=isum+icalc(i)                                             SEXI180
      if (isum.eq.0) goto 9003                                          SEXI181
                                                                        SEXI182
c --- calculation of shell if icalc(ia2) = 1 --                         SEXI183
      do 0002 ia2=1,iatom                                               SEXI184
       if (icalc(ia2).eq.1) then                                        SEXI185
c ---    cluster expansion begins --                                    SEXI186
         do 0001 i=ipr,ipr+1                                            SEXI187
            write(i,'(/a/)')'  Starting cluster expansion'              SEXI188
 0001    continue                                                       SEXI189
                                                                        SEXI190
c ---    expand the unit cell cluster around central atom --            SEXI191
         call expand (ia2,iatom,ipt,st,icntra(ia2),cel,cosi,            SEXI192
     &                itot,dmax,did,iortho,t)                           SEXI193
c ---    heap sort atom distances --                                    SEXI194
         call heapst (itot,did,4)                                       SEXI195
c ---    group into shells and write shell table ---                    SEXI196
         ibook=0                                                        SEXI197
         call shelgr (condis,conang,itot,did,elemt,ide,isearc,          SEXI198
     &                ia2,ish,dist,nshell,isub,ibook,icomp)             SEXI199
c ---    print out the vertices list for each shell                     SEXI200
         call priver(nshell,ish,ia2,elemt,ide,did,isub,icomp,title,igr, SEXI201
     &               ispa)                                              SEXI202
       end if                                                           SEXI203
 0002 continue                                                          SEXI204
      goto 9999                                                         SEXI205
                                                                        SEXI206
c --- some error labels --                                              SEXI207
 9001 continue                                                          SEXI208
      write (6,'(/a)')'  No atom card found - calculation finished '    SEXI209
      goto 9999                                                         SEXI210
 9002 continue                                                          SEXI211
      write (*,'(a,a,f7.4)') '  Cluster expansion omitted - d(max)',    SEXI212
     &                       ' too low -',dmax                          SEXI213
      goto 9999                                                         SEXI214
 9003 continue                                                          SEXI215
      write (6,'(/a)')'  No shells flagged to calculate -  finished '   SEXI216
      goto 9999                                                         SEXI217
                                                                        SEXI218
c --- all done -                                                        SEXI219
 9999 continue                                                          SEXI220
      write(*,'(/a/)') '  Bye bye from SEXIE - have a nice day'         SEXI221
      close (unit=ipr+1,status='keep')                                  SEXI222
      end                                                               SEXI223
