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