      subroutine geomet (nshell,did,ipoly,conang,icr,                   GEOM001
     &                   iclust,idim,ish,ims,averag,disoct)             GEOM002
c ----------------------------------------------------------------------GEOM003
c     This program module is part of the author's X-ray Library         GEOM004
c     and may not be sold or reproduced for commercial purposes.        GEOM005
c     Although it has been tested carefully, there is no warranty       GEOM006
c     regarding accuracy or freedom from side-effects.                  GEOM007
c         (C) 1990 Bernhard Rupp and Bryan Smith                        GEOM008
c         Major revision 4.1 : B.Rupp                                   GEOM009
c ----------------------------------------------------------------------GEOM010
c     GEOMET finds possible coordination geometry (polyhedron type)     GEOM011
c ----------------------------------------------------------------------GEOM012
      parameter (iamax=40,ispa=6)                                       GEOM013
                                                                        GEOM014
      dimension did(iclust,idim),xo(3),x1(3),x2(3),an(iamax),nan(ispa)  GEOM015
      dimension spa(ispa),ish(ims)                                      GEOM016
                                                                        GEOM017
c --- spa(ispa) contains special angle data --                          GEOM018
      data (spa(i),i=1,ispa)/                                           GEOM019
     & 109.4712, 90.0, 180.0, 60.0, 120.0, 70.5288/                     GEOM020
                                                                        GEOM021
      ishell=ish(nshell)                                                GEOM022
                                                                        GEOM023
c --- retrieve central atom cart.coord --                               GEOM024
      do 0001 i=1,3                                                     GEOM025
 0001    xo(i)=did(1,i+5)                                               GEOM026
                                                                        GEOM027
c --- must set ibeg to 1 because first loop skips --                    GEOM028
      ibeg=1                                                            GEOM029
c --- calculate where shell begins in did --                            GEOM030
      do 0101 i=1,nshell-1                                              GEOM031
         ibeg=ibeg+ish(i)                                               GEOM032
 0101 continue                                                          GEOM033
      ibeg=ibeg+1                                                       GEOM034
                                                                        GEOM035
c --- permutate points for angles --                                    GEOM036
      an(1)=0.0                                                         GEOM037
      iang=0                                                            GEOM038
c --- ishell-1 permutations --                                          GEOM039
      iper=ishell-1                                                     GEOM040
      do 0002 j=0,iper                                                  GEOM041
         do 0003, k=j+1,iper                                            GEOM042
            if (iang.ge.iamax) then                                     GEOM043
c ---       DON'T even TRY to find coordination geometry ---            GEOM044
                ipoly=icr-1                                             GEOM045
                return                                                  GEOM046
            end if                                                      GEOM047
            do 0004 l=1,3                                               GEOM048
               x1(l)=did(ibeg+j,l+5)                                    GEOM049
               x2(l)=did(ibeg+k,l+5)                                    GEOM050
 0004       continue                                                    GEOM051
            call vecang(xo,x1,xo,x2,ang)                                GEOM052
            do 1003 i1=1,iang                                           GEOM053
               if (abs(ang-an(i1)).lt.conang) goto 2003                 GEOM054
 1003       continue                                                    GEOM055
c ---       else record NEW angle ---------                             GEOM056
            iang=iang+1                                                 GEOM057
            an(iang)=ang                                                GEOM058
 2003       continue                                                    GEOM059
 0003    continue                                                       GEOM060
 0002 continue                                                          GEOM061
                                                                        GEOM062
c --- search and flag (nan) the set of special angles (spa) ---         GEOM063
      do 3000 i=1,ispa                                                  GEOM064
         nan(i)=0                                                       GEOM065
 3000 continue                                                          GEOM066
      iss=0                                                             GEOM067
c --- check for special angles --                                       GEOM068
      do 3001 i=1,iang                                                  GEOM069
c --- take all angles --                                                GEOM070
         do 3001 j=1,ispa                                               GEOM071
c ---    check em against all special angles ---                        GEOM072
            if (abs(an(i)-spa(j)).lt.conang) then                       GEOM073
c ---          mark the special angle j as found if first occurence --  GEOM074
               if (nan(j).eq.0) then                                    GEOM075
                   nan(j)=1                                             GEOM076
c ---              increase the special angle counter ---               GEOM077
                   iss=iss+j                                            GEOM078
               end if                                                   GEOM079
            endif                                                       GEOM080
 3001 continue                                                          GEOM081
                                                                        GEOM082
c --- planar ? -                                                        GEOM083
      if (ishell.gt.2) then                                             GEOM084
         call planar (did,iplan,ibeg,ishell,conang,iclust,idim,averag)  GEOM085
      endif                                                             GEOM086
                                                                        GEOM087
c --- possible arrangements for each coord.number                       GEOM088
      if (ishell.eq.1) then                                             GEOM089
         ipoly=1                                                        GEOM090
c        - one other atom point -                                       GEOM091
      elseif (ishell.eq.2) then                                         GEOM092
         if (iss.eq.3) then                                             GEOM093
            ipoly=2                                                     GEOM094
c           - linear -                                                  GEOM095
         else                                                           GEOM096
            ipoly=3                                                     GEOM097
c           - two other atoms -                                         GEOM098
         end if                                                         GEOM099
      elseif (ishell.eq.3) then                                         GEOM100
         if (iplan.eq.1) then                                           GEOM101
            ipoly=4                                                     GEOM102
c           - planar triangle -                                         GEOM103
            if (iss.eq.5) ipoly=5                                       GEOM104
c           - regular triangle -                                        GEOM105
         else                                                           GEOM106
           ipoly=6                                                      GEOM107
c          - any three other points                                     GEOM108
         end if                                                         GEOM109
      elseif (ishell.eq.4) then                                         GEOM110
         if ((iss.eq.1).and.(iang.eq.1)) then                           GEOM111
            ipoly=7                                                     GEOM112
c           - tetrahedron -                                             GEOM113
         elseif(iss.eq.5.and.iplan.eq.1) then                           GEOM114
            ipoly=8                                                     GEOM115
c           - quadratic planar                                          GEOM116
         elseif(iss.eq.3.and.iang.eq.3) then                            GEOM117
            ipoly=9                                                     GEOM118
c           - planar rectangular -                                      GEOM119
         elseif (iplan.eq.1) then                                       GEOM120
            ipoly=10                                                    GEOM121
c           - planar fourfold -                                         GEOM122
         else                                                           GEOM123
            ipoly=11                                                    GEOM124
c           - general fourfold -                                        GEOM125
         end if                                                         GEOM126
         if(ipoly.eq.11) then                                           GEOM127
c ---       search for (very) distorted tetrahedron ---                 GEOM128
            iss=0                                                       GEOM129
c ---       check for special angle 109 --                              GEOM130
            do 3201 i=1,iang                                            GEOM131
c ---          take all angles check em against special angle 1 ---     GEOM132
               if (abs(an(i)-spa(1)).gt.conang) iss=1                   GEOM133
 3201       continue                                                    GEOM134
            if(iss.eq.0) ipoly=39                                       GEOM135
         end if                                                         GEOM136
      elseif (ishell.eq.5) then                                         GEOM137
         if (iplan.eq.1) then                                           GEOM138
            ipoly=12                                                    GEOM139
c           - planar fivefold -                                         GEOM140
         else                                                           GEOM141
            ipoly=13                                                    GEOM142
c           - general fivefold -                                        GEOM143
         end if                                                         GEOM144
      elseif (ishell.eq.6) then                                         GEOM145
         if (iplan.eq.1) then                                           GEOM146
            if (iss.eq.12) then                                         GEOM147
               ipoly=14                                                 GEOM148
c              - regular hexagonal -                                    GEOM149
            else                                                        GEOM150
               ipoly=15                                                 GEOM151
c              - planar hexagonal -                                     GEOM152
            end if                                                      GEOM153
         elseif (iss.eq.5.and.iang.eq.2) then                           GEOM154
            ipoly=16                                                    GEOM155
c           - octahedral -                                              GEOM156
         elseif (nan(3).eq.0.and.iang.eq.3) then                        GEOM157
            ipoly=18                                                    GEOM158
c           - trigonal prismatic -                                      GEOM159
         elseif (nan(3).eq.1.and.iang.eq.3) then                        GEOM160
            ipoly=19                                                    GEOM161
c           - trigonal antiprismatic -                                  GEOM162
         else                                                           GEOM163
            ipoly=20                                                    GEOM164
c           - general sixfold -                                         GEOM165
         end if                                                         GEOM166
         if (ipoly.gt.16) then                                          GEOM167
c ---       special treatment for distorted octahedron --               GEOM168
            do 3100 i=1,ispa                                            GEOM169
               nan(i)=0                                                 GEOM170
 3100       continue                                                    GEOM171
            iss=0                                                       GEOM172
c ---       this search might change the number of angles ---           GEOM173
            iiang=0                                                     GEOM174
c ---       check for special angles --                                 GEOM175
            do 3101 i=1,iang                                            GEOM176
c ---          take all angles --                                       GEOM177
               do 3101 j=1,ispa                                         GEOM178
c ---             check em against all special angles ---               GEOM179
                  if (abs(an(i)-spa(j)).lt.disoct) then                 GEOM180
c ---               mark the special angle j as found --                GEOM181
                    if (nan(j).eq.0) then                               GEOM182
                        nan(j)=1                                        GEOM183
c ---                   increase the special angle counter ---          GEOM184
                        iss=iss+j                                       GEOM185
c ---                   increase the angle counter ---                  GEOM186
                        iiang=iiang+1                                   GEOM187
                    end if                                              GEOM188
                  endif                                                 GEOM189
 3101       continue                                                    GEOM190
            iang=iiang                                                  GEOM191
            if (iss.eq.5.and.iang.eq.2) then                            GEOM192
               ipoly=17                                                 GEOM193
c           -  distorted octahedral -                                   GEOM194
            end if                                                      GEOM195
         end if                                                         GEOM196
      elseif (ishell.eq.8) then                                         GEOM197
         if (iplan.eq.1) then                                           GEOM198
            ipoly=21                                                    GEOM199
c           - planar eightfold -                                        GEOM200
         elseif (iss.eq.10.and.iang.eq.3) then                          GEOM201
            ipoly=22                                                    GEOM202
c           - cubic -                                                   GEOM203
         elseif (nan(3).eq.1.and.iang.eq.5) then                        GEOM204
            ipoly=23                                                    GEOM205
c           - quadratic prismatic -                                     GEOM206
         elseif (nan(3).eq.0.and.iang.eq.4) then                        GEOM207
            ipoly=24                                                    GEOM208
c           - quadratic antiprismatic -                                 GEOM209
         elseif (iang.lt.16) then                                       GEOM210
            ipoly=25                                                    GEOM211
c           - prismatic -                                               GEOM212
         elseif (nan(3).eq.1.and.iang.eq.4) then                        GEOM213
            ipoly=26                                                    GEOM214
c           - antiprismatic -                                           GEOM215
         else                                                           GEOM216
            ipoly=27                                                    GEOM217
c           - general eightfold -                                       GEOM218
         end if                                                         GEOM219
      elseif (ishell.eq.12) then                                        GEOM220
         if (iss.eq.14.and.iang.eq.4) then                              GEOM221
            ipoly=28                                                    GEOM222
c           - simple cubeoctahedron -                                   GEOM223
         elseif (nan(3).eq.1.and.iang.eq.7) then                        GEOM224
            ipoly=29                                                    GEOM225
c           - regular hexagonal prism -                                 GEOM226
         elseif (nan(3).eq.0.and.nan(1).eq.0.and.iang.eq.9) then        GEOM227
            ipoly=30                                                    GEOM228
c           - irregular hexagonal prism -                               GEOM229
         elseif (nan(3).eq.0.and.nan(1).eq.0.and.iang.eq.4) then        GEOM230
            ipoly=31                                                    GEOM231
c           - truncated tetrahedron -                                   GEOM232
         elseif (nan(3).eq.0.and.nan(1).eq.0.and.iang.eq.5) then        GEOM233
            ipoly=32                                                    GEOM234
c           - edge and corner truncated tetrahedron -                   GEOM235
         else                                                           GEOM236
            ipoly=icr-1                                                 GEOM237
         end if                                                         GEOM238
      elseif (ishell.eq.24) then                                        GEOM239
         if (nan(3).eq.1.and.nan(1).eq.1.and.iang.eq.11) then           GEOM240
            ipoly=34                                                    GEOM241
c           - edge truncated cube -                                     GEOM242
         elseif (nan(3).eq.1.and.iang.lt.13) then                       GEOM243
            ipoly=33                                                    GEOM244
c           - hexa-cubeoctahedron -                                     GEOM245
         else                                                           GEOM246
            ipoly=icr-1                                                 GEOM247
         end if                                                         GEOM248
      elseif (ishell.eq.30) then                                        GEOM249
         if (nan(1).eq.1.and.nan(2).eq.1) then                          GEOM250
            ipoly=35                                                    GEOM251
c           - truncated cube w/ prodruding octahedron -                 GEOM252
         else                                                           GEOM253
            ipoly=icr-1                                                 GEOM254
         end if                                                         GEOM255
      elseif (ishell.eq.32) then                                        GEOM256
         if (nan(1).eq.1.and.nan(2).eq.0) then                          GEOM257
            ipoly=36                                                    GEOM258
c           - edge truncated cube with inscribed cube -                 GEOM259
         else                                                           GEOM260
            ipoly=icr-1                                                 GEOM261
         end if                                                         GEOM262
      elseif (ishell.eq.48) then                                        GEOM263
         if (nan(3).eq.1.and.iang.le.29) then                           GEOM264
            ipoly=37                                                    GEOM265
c           - edge and corner truncated cube -                          GEOM266
         else                                                           GEOM267
            ipoly=icr-1                                                 GEOM268
         end if                                                         GEOM269
      else                                                              GEOM270
            ipoly=icr-1                                                 GEOM271
      end if                                                            GEOM272
                                                                        GEOM273
c --- check for larger plane ---                                        GEOM274
      if (iplan.eq.1.and.ishell.gt.8) then                              GEOM275
         ipoly=38                                                       GEOM276
c        - plane -                                                      GEOM277
      endif                                                             GEOM278
      return                                                            GEOM279
      end                                                               GEOM280
