1      subroutine sym_map(itype, tags, symops, nops, nareal, map, xnew,
2     &     acord, oprint, ilbl, threquiv)
3C$Id$
4      implicit none
5#include "errquit.fh"
6c
7      integer itype
8      integer nareal
9      character*16 tags(nareal)
10      integer nops
11      integer maxops
12      double precision tol
13      double precision threquiv
14      parameter(maxops=192)
15      double precision symops(maxops*3,4),xnew(nops,3),acord(3,nareal)
16c
17      integer map(nops,nareal)
18      integer ilbl(nareal)
19      logical oprint
20      integer i, j, k, ii, iop, iopnum, num
21      double precision sum, dx, dy, dz
22c
23c     creates an equiv. table for mapping of atoms under sym ops
24c     stores the matrix representations of the group operators, number
25c     of operators and equiv.
26c
27c     table format:
28c     atm #1    atm #2  ..........   atm #N
29c     op#1
30c
31c     op#2
32c     .
33c     .
34c     .
35c
36c     op#N
37c
38c     therefore, searching down a col. will give you all the atoms equivalent
39c     to atom N under all operations of the group.
40c
41c     ach 10/14/94
42c
43c     Initialize map to invalid center number so that can check
44c     all entries are set at the end
45c
46      do i = 1, nareal
47         do iop = 1, nops
48            map(iop,i) = 0
49         enddo
50      enddo
51c
52c     set tolerance
53c
54      tol=threquiv
55c
56c     symops*coordinates
57c
58      sum=0.0d+00
59      do 300 i=1,nareal
60         iopnum=0
61         do 305 iop=1,nops
62            do 310 j=1,3
63               iopnum=iopnum+1
64               do 320 k=1,3
65                  sum=sum+symops(iopnum,k)*acord(k,i)
66 320           continue
67               xnew(iop,j)=sum+symops(iopnum,4)
68               sum=0.0d+00
69 310        continue
70c
71c--   > shift it into the home cell. This operation depends on sysytem type
72c--   > in place for 3-d and molecules now.
73c
74            if(itype.eq.3) then
75 3000          do 380 ii=1,3
76                  if(dabs(xnew(iop,ii)).lt.1.0d-10) then
77                     xnew(iop,ii)=0.0d0
78                  endif
79                  if(dabs((dabs(xnew(iop,ii))-1.0d0)).lt.1.0d-10) then
80                     xnew(iop,ii)=0.0d0
81                  endif
82                  if (xnew(iop,ii).lt.(0.0d0)) then
83                      xnew(iop,ii)=xnew(iop,ii)+1.0d0
84                      goto 3000
85                  elseif (xnew(iop,ii).ge.1.0d0) then
86                     xnew(iop,ii)=xnew(iop,ii)-1.0d0
87                     goto 3000
88                  endif
89 380           continue
90            elseif(itype.eq.2) then ! surfaces
91 3020         do 382 ii=1,2
92                if(dabs(xnew(iop,ii)).lt.1.0d-10) then
93                  xnew(iop,ii)=0.0d0
94                endif
95                if(dabs((dabs(xnew(iop,ii))-1.0d0)).lt.1.0d-10) then
96                   xnew(iop,ii)=0.0d0
97                endif
98                if(xnew(iop,ii).lt.0.0d0) then
99                  xnew(iop,ii)=xnew(iop,ii)+1.0d0
100                  goto 3020
101                elseif (xnew(iop,ii).ge.1.0d0) then
102                  xnew(iop,ii)=xnew(iop,ii)-1.0d0
103                  goto 3020
104                endif
105 382         continue
106            endif
107 305     continue
108c--   > load tags equiv. to the tag of atom i for search
109         num=0
110         do 400 j=1,nareal
111            if(tags(i).eq.tags(j)) then
112               num=num+1
113               ilbl(num)=j
114            endif
115 400     continue
116c--   > search those tags to identify atoms in xnew
117         do 410 j=1,nops
118            do 420 k=1,num
119               dx=dabs(acord(1,ilbl(k))-xnew(j,1))
120               if(dx.lt.tol) then
121                  dy=dabs(acord(2,ilbl(k))-xnew(j,2))
122                  if(dy.lt.tol) then
123                     dz=dabs(acord(3,ilbl(k))-xnew(j,3))
124                     if(dz.lt.tol) then
125                        map(j,i)=ilbl(k)
126                        goto 421
127                     endif
128                  endif
129               endif
130 420        continue
131            write(6,*) ' atom = ', i
132            write(6,*) ' op   = ', j
133            write(6,*) ' xnew   ', (xnew(j,k),k=1,3)
134            call util_flush(6)
135            call errquit('sym_map: no match for 1000*op+atom ',1000*j+i,
136     &       INPUT_ERR)
137 421        continue
138 410     continue
139 300  continue
140c
141c     Check the map has all entries in it.
142c
143      do i = 1, nareal
144         do j = 1, nops
145            if (map(j,i).le.0 .or. map(j,i).gt.nareal) call errquit
146     $           ('sym_map: invalid map for 1000*op+atom ',1000*j+i,
147     &       INPUT_ERR)
148         enddo
149      enddo
150c
151c--   > test prints (prints map)
152c
153c     write(*,*) 'map printed out nops rows nareal columns'
154c     do 500 i=1,nops
155c     write(*,'(8(i2,3x))') (map(i,j),j=1,nareal)
156c     500 continue
157c
158      if (oprint) write(*,10)
159 10   format(//,18x,'<<< SYMMETRY OPERATION/ATOM MAPPING BUILT >>>')
160      return
161      end
162
163
164
165
166