                          PROGRAM SEXIE                                 SEXI001
	use portlib
c ----------------------------------------------------------------------SEXI002
c     This program module is part of the author's 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 FORTRAN version 5.0b                     SEXI009
c ----------------------------------------------------------------------SEXI010
c  A. Input file : SEXIE.INP (unit=5)                                   SEXI011
c                                                                       SEXI012
c    LINE 1: A(60)                                                      SEXI013
c     title  : title line , 60 characters                               SEXI014
c    LINE 2: F,F,I,F,I                                                  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 search (1=on)                  SEXI019
c     conang : range within angles are considered to be special angles  SEXI020
c     disoct : cutoff angle for recognition of distorted octahedra      SEXI021
c     iwater : if gt.0 optinal line 2A has to follow                    SEXI022
c    LINE 2A : iwater=1                                                 SEXI023
c              izell : flag for full unit cell output                   SEXI024
c              lo,hi,loang,hiang : limits for H-bond distance and anglesSEXI025
c              omit  : do not search and print beyond first shell       SEXI026
c              bump  : bad contacts below                               SEXI027
c              ipakin: Sean Parkins weired notation for residues appliesSEXI028
c    LINE 2A : iwater=2                                                 SEXI029
c              parameters of principal axis transformation for local    SEXI030
c              charge symmetry determination                            SEXI031
c              ipower (crystal field term, 0..9)                        SEXI032
c              icharg (=1 if charge to be used)                         SEXI033
c    LINE 2A : iwater=3                                                 SEXI034
c              structure factor calculation to follow                   SEXI035
c              iano   : set to 1 for anomalous dispersion correction    SEXI036
c              lambda : guesswhat                                       SEXI037
c    LINE 3: 6(F)                                                       SEXI038
c     cel(1-6) : a,b,c,al,be,ga : cell constants                        SEXI039
c    LINE 4: A(10,I4)                                                   SEXI040
c     spcgrp : space group symbol and number (see INSHELL for details)  SEXI041
c    LINES 5 to EOF : fixed (A2,2x,A4,5f8.5,I3)                         SEXI042
c     elemt  : atom name (iat)                                          SEXI043
c     ide    : atom label (iat)                                         SEXI044
c     x,y,z  : atom positions (iat)                                     SEXI045
c     foc    : occupancy factor (iat)                                   SEXI046
c     bf     : B-factor (iwater=3) or charge (iwater=2)                 SEXI047
c     icalc  : flag for calculation of shells for corresp. atom (iat)   SEXI048
c              or ordinal number of element if sfac calc.               SEXI049
c                                                                       SEXI050
c    COMMENTS: each line preceded by a asterix * is treated as a        SEXI051
c              comment and ignored.                                     SEXI052
c                                                                       SEXI053
c  A1: Input file SEXIE.HKL : free format generated indices for sfac    SEXI054
c              calculation : h,k,l,d,sin(L)/d                           SEXI055
c                                                                       SEXI056
c                                                                       SEXI057
c  B. Variables used in program :                                       SEXI058
c     iatom  : number of inequivalent atom positions  (max iat)         SEXI059
c     ns     : number of point transformations (max 192)                SEXI060
c     ts     : translation vector of equipoint transformation           SEXI061
c     rs     : rotation martix of equipoint transformation              SEXI062
c     imult  : multiplicity of position of atom(iat)                    SEXI063
c     igen   : multiplicity of general position                         SEXI064
c     iortho : flag for orthogonal systems 1                            SEXI065
c     iwarn  : contains warning flag if equipoint atoms coincide        SEXI066
c     ucell  : array stores all generated atom positions in unit cell   SEXI067
c              its kind and and the number of its symmetry operator     SEXI068
c     iatcnt : number of atoms in ucell                                 SEXI069
c     rofats : number of each kind of atoms in ucell corrected for focc SEXI070
c     did    : position, distance, atom code, cart.coord, symmetry op.  SEXI071
c              of atoms around a central atom within dmax               SEXI072
c     alatom : contains all atoms in cluster around unit cell centre    SEXI073
c     itotat : is the total number of atoms in array alatom             SEXI074
c     ialsym : array containing the symmetry operators of alatom atoms  SEXI075
c     cosi,sinu,cosq : cosine, sine and cosine**2 of cell angles        SEXI076
c     t      : coordinate transformation matrix into cartesian coord.   SEXI077
c     ihisto : distance histogram, also histo,ibhist                    SEXI078
c                                                                       SEXI079
c  C. Variables used in mkshel and priver only :                        SEXI080
c     ishell : coord.number of shell(ish)                               SEXI081
c     nshell : total number of coord. shells within dmax                SEXI082
c     ish    : array containing number of atoms (ishell) for each shell SEXI083
c     iss2   : contains the atom type for each subshell(is)             SEXI084
c     igeom  : stores geometry for each shell(ish)                      SEXI085
c                                                                       SEXI086
c  D. Variables needed for printing :                                   SEXI087
c     icomp  : type of computer compiled for                            SEXI088
c     lv     : line vertical character for boxes                        SEXI089
c     an     : Angstroem character                                      SEXI090
c     de     : degree character                                         SEXI091
c     pm     : plus/minus sign character                                SEXI092
c     hf     : 1/2 character                                            SEXI093
c ----------------------------------------------------------------------SEXI094
c     This program compromises between memory usage and computation     SEXI095
c     speed so that it can run on 486 PC's as well. Please take care    SEXI096
c     when changing assignment of variables which might be indexing     SEXI097
c     variables (arrays) connected to primary data arrays.              SEXI098
c     For easier modifications, no common blocks are used and all       SEXI099
c     global variables are explicitly referenced by name or value       SEXI100
c ----------------------------------------------------------------------SEXI101
c     Any improvements and critizism (in particular for a fast          SEXI102
c     mathematical algorithm to find the coordination geometry)         SEXI103
c     should be sent to :                                               SEXI104
c     Bernhard Rupp                                                     SEXI105
c     RUPP1@LLNL.GOV                                                    SEXI106
c     LLNL - BBRP                                                       SEXI107
c ----------------------------------------------------------------------SEXI108
c     Parameters used in main program and subroutines:                  SEXI109
c     iat    : max. number of nonequivalent sites or atoms in cell      SEXI110
c     iclust : nax. number of atoms in cluster                          SEXI111
c     istor  : max. number of atoms in unit cell                        SEXI112
c     ims    : max. number of shells (100)                              SEXI113
c     ipr    : print file unit 6 - do not change                        SEXI114
c     idim   : number of elements in did (leave at 9)                   SEXI115
c     ialdim : 3 coordinates plus atom type and generator in alatom     SEXI116
c     igeo   : 41 recognized geometries in geom                         SEXI117
c     incel  : 3 coordinates plus type plus operator in ucell           SEXI118
c     nel    : number of recognized chemical element symbols (105)      SEXI119
c ----------------------------------------------------------------------SEXI120
c     files opened : 3 = infeff.xyz, coordinat file for FEFF & LOCSYM   SEXI121
c                    4 = configuration file SEXIE.CFG                   SEXI122
c                    5 = sexie.inp                                      SEXI123
c                    6 = screen                                         SEXI124
c                    7 = sexie.out, standard output file                SEXI125
c                    8 = xyz.dat and atom.dat, contain symmtry operatorsSEXI126
c                    8 = water.out in main routine                      SEXI127
c                    9 = cell.dat, contains unit cell data via CELOUT   SEXI128
c                   10 = jobnam.hkl, index file for sfac calculation    SEXI129
c                   17 = cormat.dat, coordinate orthogonalization matrixSEXI130
c                   20 = dhistogr.dat, radial dist.distribution functionSEXI131
c ----------------------------------------------------------------------SEXI132
c --- generic data parameters --                                        SEXI133
      parameter (iat=2000,istor=20000,ims=100,ipr=6,iclust=20000,idim=9)SEXI134
      parameter (ialdim=4,igeo=41,incel=5,nel=105)                      SEXI135
	integer*4 time_zero
c --- these for histograms --                                           SEXI136
      parameter (nhisto=78)                                             SEXI137
      character elemt*2,ide*4,title*60,revdat*8,icomp*4,avern*3         SEXI138
      character lv,an,de,pm,hf,geom*25                                  SEXI139
c --- water polygon search --                                           SEXI140
      parameter (iri=100)                                               SEXI141
      integer irff(iri),itetra(iat),ihyd(20,2),ihydtp(iat,2)            SEXI142
      Integer iwarr(iat,2)                                              SEXI143
      real hi,lo,hiang,loang,rindid(iri*7,6),ringsm(iri),watat(iat,4)   SEXI144
c --- bad contacts --                                                   SEXI145
      parameter (ibadct=50)                                             SEXI146
      real omit, bad(ibadct,3)                                          SEXI147
c --- atom input related --                                             SEXI148
      dimension elemt(iat),ide(iat),x(iat),y(iat),z(iat),igr(36)        SEXI149
      dimension bf(iat),foc(iat),iwarn(iat),icalc(iat),imult(iat),cel(6)SEXI150
c --- space group related ---                                           SEXI151
      dimension ts(3,24),fs(3,3,24),cosi(3),sinu(3),t(3,3),cosq(3)      SEXI152
c --- mass storage....                                                  SEXI153
      integer*2 ialsym(iclust*10,4)                                     SEXI154
      dimension did(iclust,idim),ucell(istor,incel)                     SEXI155
      dimension alatom(iclust*10,ialdim)                                SEXI156
c --- intermediate sorting storage ---                                  SEXI157
      dimension rheap(idim)                                             SEXI158
c --- shell grouping business                                           SEXI159
      dimension ish(ims),igeom(ims),geom(igeo)                          SEXI160
c --- the following for radial distribution function histogram --       SEXI161
      dimension ihisto(nhisto),histo(nhisto),ibhist(451)                SEXI162
c --- structure factor calculation related ---                          SEXI163
      character jobnam*60                                               SEXI164
      real rofats(nel),f000,lambda,xdens                                SEXI165
                                                                        SEXI166
                                                                        SEXI167
c --- block data for                                                    SEXI168
      data (geom(i),i=1,igeo)/                                          SEXI169
     & 'one other atom','linear','general twofold',                     SEXI170
     & 'planar triangle','regular triangle','general trifold',          SEXI171
     & 'tetrahedral','quadratic planar','plane rectangle',              SEXI172
     & 'planar fourfold','general fourfold',                            SEXI173
     & 'planar fivefold','general fivefold',                            SEXI174
     & 'regular hexagon','planar sixfold','octahedral',                 SEXI175
     & 'distorted octahedral','trigonal prism',                         SEXI176
     & 'trigonal antiprism','general sixfold','planar eightfold',       SEXI177
     & 'cubic','quadratic prism','quadratic antiprism',                 SEXI178
     & 'general prism','general antiprism','general eightfold',         SEXI179
     & 'simple cuboctahedron','reg. hexagonal prism','irreg. hex. prism'SEXI180
     & ,'truncated tetrahedron','doubly trnctd tetrahedron',            SEXI181
     & 'hexa-cuboctahedron', 'edge truncated cube',                     SEXI182
     & 'trnctd cube + octahedron','edge trnctd cube + cube',            SEXI183
     & 'doubly truncated cube','planar','distorted tetrahedron',        SEXI184
     & 'no special name','    -N/A- '/                                  SEXI185
                                                                        SEXI186
c --- revision data -- PUT ANY REVISION MARKS IN HERE PLEASE ---------  SEXI187
      revdat='08/21/95'                                                 SEXI188
      avern='5.0'
	time_zero=time()                                                  SEXI189
                                                                        SEXI190
c --- REVISION 5.0 :                                                    SEXI191
c     input file contains now foc, bf (see above for format details)    SEXI192
c     ATOM LINES (line 5 to EOF) : fixed (A2,2x,A4,5f8.5,I3)            SEXI193
c     call to GETINP formal parameter bf(iat) and foc(iat) added        SEXI194
c     anomalous structure factor calculation added                      SEXI195
c     infeff.xyz includes focc now (xyzout.for)                         SEXI196
c                                                                       SEXI197
c --- REVISION 4.2 :                                                    SEXI198
c     principal axis transformation included for symmetry determination SEXI199
c                                                                       SEXI200
c --- REVISION 4.1 :                                                    SEXI201
c     Major changes to storage method :                                 SEXI202
c     all atoms in unit cell IUCELL now stored in UCELL                 SEXI203
c     all dimensions handed down into subroutines for easy change       SEXI204
c     number of expansions calculated beforehand                        SEXI205
c     statistics of the computation given as output                     SEXI206
c     no more search for central atom close to origin                   SEXI207
c     histogram of distances (rad. distribution function) included      SEXI208
c     search for all atoms around unit cell within dmax first (ALATOM)  SEXI209
c     leads to shorter search in each cycle (important for proteins)    SEXI210
c     dimensions can be changed now to accomodate even proteins         SEXI211
c                                                                       SEXI212
c --- REVISION 4.0 :                                                    SEXI213
c     input file contains now charge (see above for format details)     SEXI214
c     ATOM LINES (line 5 to EOF) : fixed (A4,A4,4f8.5,I3)               SEXI215
c     call to GETINP formal parameter charge (iat) added                SEXI216
c     subroutine system renamed to xtlsys (UNIX R6000 problem)          SEXI217
c     subroutine tensor prints coordinate orthogonalization matrix      SEXI218
c     into cormat.dat (unit 17)                                         SEXI219
c                                                                       SEXI220
c --- REVISION 3.0 includes Subroutine CELOUT, PRIVER :                 SEXI221
c     numerical error in transformation from fractional to cartesian    SEXI222
c     coordinates - fixed                                               SEXI223
c     basis tramsformation made standard with IUCr                      SEXI224
c     I41/amd and P3221 decoded wrong due to rounding precision-fixed   SEXI225
c     all system files upper case now                                   SEXI226
c     PRIVER.FOR reads ispa now -  A.T.M.                               SEXI227
                                                                        SEXI228
c --- initialize configuration --                                       SEXI229
      call inicfg(icomp,lv,an,de,pm,hf)                                 SEXI230
                                                                        SEXI231
c --- defauft unit screen = 6, output file = 7 -                        SEXI232
      open (unit=ipr+1,file='sexie.out',status='unknown')               SEXI233
                                                                        SEXI234
c --- reset transformation matrix, space group symbol and div. flags    SEXI235
      call reset (ts,fs,igr,isyst,ns,isymce,iwarn,icalc,iat,nhisto,     SEXI236
     &            histo)                                                SEXI237
                                                                        SEXI238
c --- read input file --                                                SEXI239
      call getinp (title,elemt,ide,x,y,z,cel,iatom,igr,isearc,iwater,   SEXI240
     &conang,condis,dmax,ispa,icalc,foc,bf,iat,disoct,iano,lambda,      SEXI241
     &izell,lo,hi,loang,hiang,omit,bump,ipakin,ifiwnr,ipower,icharg,    SEXI242
     &jobnam)                                                           SEXI243
                                                                        SEXI244
c --- reset if unrealistic parameters for shell search -                SEXI245
      if (condis.lt.0.00001) condis=0.00001                             SEXI246
      if (conang.lt.0.05) conang=0.05                                   SEXI247
                                                                        SEXI248
c --- check condis                                                      SEXI249
      if (condis.gt.dmax) condis=dmax                                   SEXI250
                                                                        SEXI251
c --- add condis to dmax                                                SEXI252
      dmax=dmax+condis                                                  SEXI253
                                                                        SEXI254
c --- retrieve crystal system from cell constants -                     SEXI255
      call xtlsys (cel,isystm,iortho)                                   SEXI256
                                                                        SEXI257
c --- calculate trigonometric terms for storage -                       SEXI258
      call trigo(cel,cosi,sinu,cosq)                                    SEXI259
                                                                        SEXI260
c --- calculate cell volume  -                                          SEXI261
      call volume (cel,cosq,cosi,vez,iortho)                            SEXI262
                                                                        SEXI263
c --- calculate the transformation matrix T for cell to orthonormal -   SEXI264
      call trigo2 (cosi,sinu,t,cosq)                                    SEXI265
                                                                        SEXI266
c --- decode space group for equipoint transformations                  SEXI267
      call equipt (isyst,isymce,isymbr,ns,ts,fs,igr)                    SEXI268
                                                                        SEXI269
c --- ns = number of point symmetry operations --                       SEXI270
                                                                        SEXI271
c --- get number of bravais lattice -                                   SEXI272
      ibravl=ibrav(isymbr)                                              SEXI273
                                                                        SEXI274
c --- print first part of output while checking consistency of system   SEXI275
      call oneout (title,dmax,conang,condis,isearc,iatom,vez,cel,       SEXI276
     &isystm,isyst,igr,isymce,ibravl,ispa,t,cosi,revdat,icomp,lv,an,de, SEXI277
     &avern,izell,lo,hi,loang,hiang,omit,iwater,bump,iano,lambda)       SEXI278
                                                                        SEXI279
c --- print equipoint matrix and decoded equipoints -                   SEXI280
      call eqiout (ts,fs,ns,isymce,ibravl,igen,icomp,lv,hf)             SEXI281
                                                                        SEXI282
c --- end program if no atom card found --                              SEXI283
      if (iatom.eq.0) goto 9001                                         SEXI284
                                                                        SEXI285
      write(*,'(/a/)')'  Expanding asymmetric unit into unit cell......'SEXI286
                                                                        SEXI287
c --- multiplicities of atom positions and print atoms in asym.cell --  SEXI288
      call multi2 (iatom,ibravl,x,y,z,fs,ts,imult,isymce,ns,igen,       SEXI289
     &             iwarn,cel,ucell,incel,istor,iat,iatcnt)              SEXI290
c      print*,'iatcnt=',iatcnt                                          SEXI291
c      pause                                                            SEXI292
                                                                        SEXI293
c --- expand the matrix T with the cell dimensions --                   SEXI294
      call trigo3(cel,t)                                                SEXI295
                                                                        SEXI296
c --- print atom positions and unit cell information -                  SEXI297
      call celout (iatom,elemt,ide,x,y,z,imult,iwarn,foc,bf,            SEXI298
     &     vez,t,cel,icomp,title,ucell,incel,istor,iat,lv,an,           SEXI299
     &     f000,rofats,nel,xdens,izell)                                 SEXI300
                                                                        SEXI301
                                                                        SEXI302
c --- end program here, if dmax lt 1.0 --                               SEXI303
      if (dmax.lt.1.0) goto 9002                                        SEXI304
                                                                        SEXI305
c --- structure factor calculation --                                   SEXI306
      if (iwater.eq.3) then                                             SEXI307
                                                                        SEXI308
c         write(*,'(/a/)')'  Begin of structure factor calculation.....'SEXI309
         call sfcalc (ucell,istor,incel,iatcnt,bf,foc,elemt,            SEXI310
     &   iat,rofats,nel,xdens,f000,lambda,iano,isymce,jobnam,icomp)     SEXI311
         goto 9004                                                      SEXI312
      end if                                                            SEXI313
c     subroutine sfcalc (ucell,istor,incel,iatcnt,bf,foc,               SEXI314
c    & elemt,iat,rofats,nel,xdens,f000,lambda,iano,isymce,jobnam,icomp) SEXI315
                                                                        SEXI316
c --- calculate last shell and end program here if no shells --         SEXI317
      isum=0                                                            SEXI318
      do 0003 i=1,iatom                                                 SEXI319
         isum=isum+icalc(i)                                             SEXI320
         if (icalc(i).gt.0) ifeff=i                                     SEXI321
 0003 continue                                                          SEXI322
      if (isum.eq.0) goto 9003                                          SEXI323
                                                                        SEXI324
c --- determine size of the unit cell cluster and find all atoms --     SEXI325
      call clusiz (icomp,cel,cosi,itotat,dmax,alatom,iortho,t,lv,iatom, SEXI326
     &     nexp,imult,iucell,ucell,incel,iclust,ialdim,istor,iat,ialsym)SEXI327
                                                                        SEXI328
c --- reset the ring found counters --                                  SEXI329
      iring=1                                                           SEXI330
      irf=0                                                             SEXI331
      irfdid=0                                                          SEXI332
      ihand=0                                                           SEXI333
      ibad=0                                                            SEXI334
      iwat=0                                                            SEXI335
c --- omit is kept in romit because it's changed for priver --          SEXI336
      romit=omit                                                        SEXI337
                                                                        SEXI338
c --- calculation of shell if icalc(ia) = 1 --                          SEXI339
      do 0002 ia=1,iatom                                                SEXI340
       if (icalc(ia).eq.1) then                                         SEXI341
c ---    cluster expansion begins --                                    SEXI342
         do 0001 i=ipr,ipr+1                                            SEXI343
            write(i,'(/a,i4,a,a2,a4/)')                                 SEXI344
     &         '  Starting to search ligands around atom Nr.',ia,       SEXI345
     &         ' : ',elemt(ia),ide(ia)                                  SEXI346
 0001    continue                                                       SEXI347
c ---    expand the unit cell cluster around central atom ia --         SEXI348
         call corshl (ia,cel,cosi,itot,dmax,did,iortho,t,alatom,        SEXI349
     &   ucell,incel,imult,istor,iat,iclust,idim,itotat,ialdim)         SEXI350
                                                                        SEXI351
c            print*,'-',ide(did(ia,5)),'-'                              SEXI352
c            if (ia.eq.167) then                                        SEXI353
c             iyes=1                                                    SEXI354
c             write(7,*)'Before heapst'                                 SEXI355
c               do 1001 ika=1,itot                                      SEXI356
c                  write(7,'(1x,a2,1x,a4,4f8.5)')                       SEXI357
c    &       elemt(did(ika,5)),ide(did(ika,5)),(did(ika,joka),joka=1,4) SEXI358
                                                                        SEXI359
c 1001          continue                                                SEXI360
c            end if                                                     SEXI361
                                                                        SEXI362
c ---    sort them according to distance (did(i,4)) --                  SEXI363
         call heapst (itot,did,4,iclust,idim,rheap)                     SEXI364
                                                                        SEXI365
c            print*,'-',ide(did(1,5)),'-'                               SEXI366
c            if (int(did(1,5)).eq.167) then                             SEXI367
c             write(7,*)'After heapst'                                  SEXI368
c               do 1002 ika=1,itot                                      SEXI369
c                  write(7,'(1x,a2,1x,a4,4f8.5)')                       SEXI370
c    &       elemt(did(ika,5)),ide(did(ika,5)),(did(ika,joka),joka=1,4) SEXI371
c                                                                       SEXI372
c 1002          continue                                                SEXI373
c            end if                                                     SEXI374
                                                                        SEXI375
                                                                        SEXI376
c       do 5000 ii=1,itot                                               SEXI377
c5000    write(*,'(a,9f8.4)')' did=',(did(ii,iii),iii=1,9)              SEXI378
c         pause                                                         SEXI379
                                                                        SEXI380
c ---    add into the histogram and check bad distances --              SEXI381
         call adhist(itot,did,iclust,idim,histo,nhisto,ibhist,          SEXI382
     &               bad,ibadct,ibad,iwater,bump)                       SEXI383
                                                                        SEXI384
c ---    group into shells and write shell table ---                    SEXI385
         ipoly=0                                                        SEXI386
         omit=romit                                                     SEXI387
         call mkshel (condis,conang,itot,did,elemt,ide,isearc,ia,       SEXI388
     &   ish,ims,nshell,icomp,iat,iclust,idim,lv,disoct,igeom,geom,igeo,SEXI389
     &   itetra,omit,ihyd,iwater)                                       SEXI390
                                                                        SEXI391
c ---    check hydration level --                                       SEXI392
         if(iwater.eq.1) then                                           SEXI393
             call hydrat(did,iclust,idim,iwat,ihydtp,iat,ihyd,itot,     SEXI394
     &       elemt,watat,ix)                                            SEXI395
             if(watat(ix,2).gt.0.001) then                              SEXI396
               omit=watat(ix,2)                                         SEXI397
             elseif(watat(ix,3).gt.0.001) then                          SEXI398
               omit=watat(ix,3)                                         SEXI399
             end if                                                     SEXI400
         end if                                                         SEXI401
                                                                        SEXI402
         omit=omit+0.001                                                SEXI403
c ---    output for infeff for last assigned shell only (bf=charge) --- SEXI404
         if (ifeff.eq.ia) then                                          SEXI405
            call xyzout (condis,itot,did,elemt,ide,ia,avern,ispa,       SEXI406
     &      title,iat,iclust,idim,bf,foc,igr)                           SEXI407
         end if                                                         SEXI408
                                                                        SEXI409
                                                                        SEXI410
c ---    print out the vertices list for each shell (bf=charge)         SEXI411
         call priver (nshell,ish,ia,elemt,ide,did,icomp,iat,iclust,     SEXI412
     &  idim,ims,lv,pm,igeom,geom,igeo,omit,ialsym,ipower,icharg,bf,foc)SEXI413
                                                                        SEXI414
c ---    do water closed loop search if iwater is set 1 ----            SEXI415
         if (iwater.eq.1) then                                          SEXI416
c            print*,'-',ide(did(1,5)),'-'                               SEXI417
            if (ide(did(1,5)).eq.'103 ') then                           SEXI418
               do 1000 ika=1,itot                                       SEXI419
                  write(7,'(1x,a2,1x,a4,4f8.5)')                        SEXI420
     &       elemt(did(ika,5)),ide(did(ika,5)),(did(ika,joka),joka=1,4) SEXI421
                                                                        SEXI422
 1000          continue                                                 SEXI423
            end if                                                      SEXI424
            call waters(itot,did,iclust,idim,iortho,cosi,iring,irf,irff,SEXI425
     &           irfdid,ihand,rindid,ringsm,iri,cel,lo,hi,loang,hiang)  SEXI426
c            call windex(itot,did,iclust,idim,iortho,cosi,iring,irf,irffSEXI427
c     &      irfdid,ihand,rindid,ringsm,iri,cel,lo,hi,loang,hiang)      SEXI428
                                                                        SEXI429
                                                                        SEXI430
         end if                                                         SEXI431
       end if                                                           SEXI432
 0002 continue                                                          SEXI433
                                                                        SEXI434
c --- extra file for the water information ---                          SEXI435
      if (iwater.eq.1) then                                             SEXI436
         open(ipr+2,file='water.out',status='unknown')                  SEXI437
                                                                        SEXI438
c ---    print ring table ---                                           SEXI439
         call rintab (ipr,irff,irf,iri,rindid,icomp,lv,geom,igeo,       SEXI440
     &   elemt,ide,iat,ihand,lo,hi,loang,hiang,itetra,condis,conang,    SEXI441
     &   bad,ibadct,ibad,bump,ihydtp,iwat,an,de,watat,iwarr)            SEXI442
      end if                                                            SEXI443
                                                                        SEXI444
c --- print the histogram                                               SEXI445
      iplot=0                                                           SEXI446
      if (iwater.eq.1) iplot=1                                          SEXI447
      call histog                                                       SEXI448
     &      (icomp,histo,ihisto,nhisto,ipr,ibhist,title,avern,iplot)    SEXI449
                                                                        SEXI450
c --- close water file --                                               SEXI451
      if (iwater.eq.1) close(ipr+2,status='keep')                       SEXI452
                                                                        SEXI453
      goto 9999                                                         SEXI454
                                                                        SEXI455
c --- some error labels --                                              SEXI456
                                                                        SEXI457
 9001 continue                                                          SEXI458
      write (6,'(/a)')'  No atom card found - calculation finished '    SEXI459
      goto 9999                                                         SEXI460
 9002 continue                                                          SEXI461
      write (*,'(a,a,f7.4)') '  Cluster expansion omitted - d(max)',    SEXI462
     &                       ' too low -',dmax                          SEXI463
      goto 9999                                                         SEXI464
 9003 continue                                                          SEXI465
      write (6,'(/a)')'  No shells flagged to calculate -  finished '   SEXI466
      goto 9999                                                         SEXI467
                                                                        SEXI468
 9004 continue                                                          SEXI469
      write (6,'(/a)')'  Structure factor calculation finished '        SEXI470
      goto 9999                                                         SEXI471
                                                                        SEXI472
c --- all done -                                                        SEXI473
 9999 continue                                                          SEXI474
      close (unit=ipr+1,status='keep') 

	write(*,'(/a,i4,a/)')' Elapsed time : ',time()-time_zero, ' sec.'

      end                                                               SEXI476
