      subroutine cluster(iatom,cel,cosi,itotat,dmax,alatom,iortho,t,    CLUS001
     &radius,nexp,ucell,incel,imult,iucell,istor,iat,iclust,ialdim,     CLUS002
     &ialsym)                                                           CLUS003
c ----------------------------------------------------------------------CLUS004
c     This program module is part of the author's X-ray Library         CLUS005
c     and may not be sold or reproduced for commercial purposes.        CLUS006
c     Although it has been tested carefully, there is no warranty       CLUS007
c     regarding accuracy or freedom from side-effects.                  CLUS008
c                   (C) 1986-1994 Bernhard Rupp                         CLUS009
c ----------------------------------------------------------------------CLUS010
c     CLUSTER generates translation vectors for cell expansion          CLUS011
c     and expands until all distances < dmax are found                  CLUS012
c ----------------------------------------------------------------------CLUS013
      integer*2 ialsym(iclust*10,4)                                     CLUS014
      dimension center(3),corna(8,3),temp(3)                            CLUS015
      dimension p(3),alatom(iclust*10,ialdim),cel(6),cosi(3)            CLUS016
      dimension xo(3),xa(3),t(3,3),ucell(istor,incel),imult(iat)        CLUS017
                                                                        CLUS018
c --- find the largest of 4 space diagonal in any oblique cell ---      CLUS019
                                                                        CLUS020
c --- create corners 000 - 111, 010 - 101, 100 - 011, and 001 - 110     CLUS021
      do 0001, i=1,3                                                    CLUS022
         center(i)=0.5                                                  CLUS023
         do 0002 j=1,8                                                  CLUS024
c ---       even                                                        CLUS025
            if(mod(j,2).eq.0) then                                      CLUS026
               corna(j,i)=1.0                                           CLUS027
            else                                                        CLUS028
               corna(j,1)=0.0                                           CLUS029
            end if                                                      CLUS030
 0002    continue                                                       CLUS031
 0001 continue                                                          CLUS032
      corna(3,2)=1.0                                                    CLUS033
      corna(4,2)=0.0                                                    CLUS034
      corna(5,1)=1.0                                                    CLUS035
      corna(6,1)=0.0                                                    CLUS036
      corna(7,3)=1.0                                                    CLUS037
      corna(8,3)=0.0                                                    CLUS038
                                                                        CLUS039
c --- make them cartesian --                                            CLUS040
      do 0004 i=1,8                                                     CLUS041
         do 0006 j=1,3                                                  CLUS042
            temp(j)=corna(i,j)                                          CLUS043
 0006    continue                                                       CLUS044
         call tracor(temp,t)                                            CLUS045
         do 0007 j=1,3                                                  CLUS046
            corna(i,j)=temp(j)                                          CLUS047
 0007    continue                                                       CLUS048
 0004 continue                                                          CLUS049
                                                                        CLUS050
      bigtit=0.0                                                        CLUS051
      do 0010 j=1,7,2                                                   CLUS052
         deltax=(corna(j+1,1) - corna(j,1))**2                          CLUS053
         deltay=(corna(j+1,2) - corna(j,2))**2                          CLUS054
         deltaz=(corna(j+1,3) - corna(j,3))**2                          CLUS055
         diag=sqrt(deltax+deltay+deltaz)                                CLUS056
         if (diag.gt.bigtit) bigtit=diag                                CLUS057
 0010 continue                                                          CLUS058
      radius=dmax+bigtit/2.0                                            CLUS059
                                                                        CLUS060
c --- find the centre of the cluster --                                 CLUS061
      do 2011 i=1,3                                                     CLUS062
         xo(i)=center(i)                                                CLUS063
 2011 continue                                                          CLUS064
                                                                        CLUS065
      itotat=0                                                          CLUS066
c --- loop over all atoms in unit cell                                  CLUS067
      do 1010 ia=1,iucell                                               CLUS068
c ---    generate the quadrant translation vectors -                    CLUS069
            do 1012 it1= -nexp,nexp                                     CLUS070
               do 1013 it2= -nexp,nexp                                  CLUS071
                  do 1014 it3= -nexp,nexp                               CLUS072
                     p(1)=it1+ucell(ia,1)                               CLUS073
                     p(2)=it2+ucell(ia,2)                               CLUS074
                     p(3)=it3+ucell(ia,3)                               CLUS075
c ---                calculate the distance to newly generated point -  CLUS076
                     do 2003 k=1,3                                      CLUS077
                         xa(k)=cel(k)*(p(k)-xo(k))                      CLUS078
 2003                continue                                           CLUS079
                     dd=distan (iortho,xa,cosi)                         CLUS080
c ---                skip if distance larger than d(max) else add --    CLUS081
                     if(dd.le.radius) then                              CLUS082
                       itotat=itotat+1                                  CLUS083
c ---                  warn if itotat exceeds iclust ---                CLUS084
                       if(itotat.gt.(iclust*10)) then                   CLUS085
                         write(*,'(a/a,i5/a)')'   ****WARNING****',     CLUS086
     &                  '  Cluster contains too many atoms :',iclust*10,CLUS087
     &                  '  You need to increase parameter ICLUST*10 '   CLUS088
                         pause                                          CLUS089
                       end if                                           CLUS090
                       do 2029 ij=1,3                                   CLUS091
                          alatom(itotat,ij)=p(ij)                       CLUS092
 2029                  continue                                         CLUS093
c ---                  compute what atom kind this is                   CLUS094
c                       inum=0                                          CLUS095
c                       do 2033 i=1,iatom                               CLUS096
c                          inum=inum+imult(i)                           CLUS097
c                          if (ia.le.inum) then                         CLUS098
c                             inum=i                                    CLUS099
c                             goto 2035                                 CLUS100
c                          end if                                       CLUS101
c 2033                  continue                                        CLUS102
c 2035                  continue                                        CLUS103
c                       id=int(ucell(ia,5))                             CLUS104
                       inum=int(ucell(ia,5))                            CLUS105
c                       if (id.ne.inum) then                            CLUS106
c                          write(*,*)                                   CLUS107
c     &                    ' ** HORROR ** problem generating equivalentsCLUS108
c                          pause                                        CLUS109
c                       end if                                          CLUS110
                       alatom(itotat,4)=inum                            CLUS111
                       ialsym(itotat,1)=it1                             CLUS112
                       ialsym(itotat,2)=it2                             CLUS113
                       ialsym(itotat,3)=it3                             CLUS114
                       ialsym(itotat,4)=int(ucell(ia,4))                CLUS115
c                       print*,(ialsym(itotat,ij),ij=1,4)               CLUS116
                     end if                                             CLUS117
 1014             continue                                              CLUS118
 1013          continue                                                 CLUS119
 1012       continue                                                    CLUS120
 1011    continue                                                       CLUS121
 1010 continue                                                          CLUS122
      return                                                            CLUS123
      end                                                               CLUS124
                                                                        CLUS125
