1cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2c
3c
4c Program gensym:
5c
6c   Given an input set of generators for any group in 3 dimensional
7c  space or less, this program will produce the matrix representations
8c  of the operators in the group. The input operators must be chosen in
9c  a particular order and according to the conventions set forth
10c  on page 729 of the International Tables of Crystallography vol. 3.
11c
12c Algorithm:
13c    The algorithm works in the following manner: given a set of
14c   generators, {G}n choose and initial generator, G(i). Multiply
15c   this generator by some operator from the set {O}m. In the first
16c   pass of the seqence O(j)=G(i) therefore O(k)=G(i)*O(j)=G(i)*G(i).
17c   Which is added to the current symmetry operation list (the matrix
18c   symops) if G(i)^2.ne.E. A new generator is then selected from the
19c   generator list and added to the bottom of symops. All previous
20c   operators in the list are then pre-multiplied by this new generator.
21c    For example, if the final generator is G3, then sympos may contain:
22c
23c   G1, G1^2, G2, G2G1, G2G1^2, G2^2, G2^2G1, G2^2G1^2, G3, G3G1,
24c   G3G1^2, G3G2, G3G2G1, G3G2G1^2, G3G2^2, G3G2^2G1, G3G2^2G1^2,
25c   G3^2, G3^2G1, G3^2G1^2, G3^2G2, G3^2G2G1, G3^2G2G1^2, G3^2G2^2,
26c   G3^2G2^2G1 and G3^2G2^2G1^2.
27c
28c    Note: Identities and repeated operators generated along the way
29c   are trapped and not added to the list, so the actual number of
30c   operators in symops may be considerably less than shown above.
31c
32c Matrix Mults:
33c    The appropriate matrix multiplications are defined in Seitz
34c   notation as:
35c
36c            {R|t}{S|u}r={R|t}(Sr+u)=RSr+Ru+t={RS}|Ru+t}r
37c
38c   R and S are normal point operations t and u are the associated
39c   translation components for whatever dimensional system you are
40c   working with. Expainations of this notation can be found in "Space
41c   Groups for Solid State Scientists", G. Burns and A.M. Glazer.
42c
43c
44c    The matrix SYMOPS contains the matrix reps. of all group operators
45c   except the identity. The variable NOPS holds the number of operators
46c   in SYMOPS.
47c
48c
49c
50c
51c                                        A.C. Hess
52c                                        D.G. Clerc
53c                                        Solid State Theory Group
54c                                        MSRC/PNL
55c                                        9/13/93
56c***********************************************************************
57      subroutine gensym(itype,numgrp,numset,symops,nops,oprint,
58     $     group_name,geom,rtdb)
59C$Id$
60      implicit real*8 (a-h,o-z)
61
62      parameter(maxops=192,tol=1.0d-07,max_gen=6)
63      character*2 kpos(-1:3),kneg(-3:1),rotoop(maxops)
64      dimension capr(3,4),caps(3,4),symops(maxops*3,4)
65*      integer indx(3)
66      double precision deter3
67      external deter3
68      dimension resop(3,4),gens(18,4),detres(3,3),cntvec(3,3)
69      character*(*) group_name
70      integer geom,rtdb
71      logical oprint,use_primitive
72      double precision s_vec(max_gen,3)
73      data kpos/' 2',' 3',' 4',' 6',' 1'/,kneg/'-1','-6','-4','-3',' m'/
74
75      logical  geom_use_primitive
76      external geom_use_primitive
77
78c
79c-->call spgen with correct system type flag to make generators
80c
81         call spgen(itype,numgrp,numset,gens,cntvec,ngen,numvec,
82     $              group_name,max_gen,s_vec,oprint)
83
84c
85c-----------------------------------------------------------------------
86c
87c--> outer loop picks up input generators (begins at 5000) inner implied
88c    loop mults. current generator by all previous operators generated
89c    and stored in symops (begins at 6000).
90c
91c--> pointers:
92c             ipos indexes the current working generator in symops
93c             igpos points at the next element to be used in the
94c                  input generator list
95c             jpos points to the current operator in the list symops
96c             nops holds the current number of operators in symops
97c
98c----------------------------------------------------------------------
99      igen=1
100      icnt1=1
101      isquare=0
102      nops=0
103      if (oprint) then
104         write(*,27)
105         write(*,29) ngen
106 27      format(/,16x,'---------------',' GROUP GENERATORS ','---------'
107     $        ,'------')
108 28      format(/,23x,'GROUP NUMBER AND NAME: ',a12)
109 29      format(/,22x,i1,' GENERATORS USED TO FORM THE GROUP')
110      endif
1115000  ipos=(icnt1-1)*3+1
112      igpos=(igen-1)*3+1
113c
114c--> pickup an input generator and store it in the symop list
115c--> set first operator equal to this generator, G(i).
116c
117      do 200 i=1,3
118         do 220 j=1,4
119            symops(ipos+(i-1),j)=gens(igpos+(i-1),j)
120            capr(i,j)=gens(igpos+(i-1),j)
121220      continue
122200   continue
123      nops=nops+1
124c
125c--> compute trace and determinant of the generator
126c
127      do 560 i=1,3
128         do 570 j=1,3
129            detres(i,j)=capr(i,j)
130570      continue
131560   continue
132      trace=0.0d+00
133      do 580 i=1,3
134         trace=trace+capr(i,i)
135580   continue
136*
137      det = deter3(detres)
138*      call ludcmp(detres,3,3,indx,det)
139*      do 590 i=1,3
140*         det=det*detres(i,i)
141*590   continue
142c
143c--> store type information
144c
145c**********************************************************************
146c  Screen out all point groups from the next calc   DGC 3/10/94
147c**********************************************************************
148      if(itype.ge.1.and.itype.le.3) then
149         itrace=idint(trace)
150c
151         if(det.lt.0.0d+00) then
152            rotoop(nops)=kneg(itrace)
153            if (oprint) write(*,30) kneg(itrace),(s_vec(igen,j),j=1,3)
154         else
155            rotoop(nops)=kpos(itrace)
156            if (oprint) write(*,30) kpos(itrace),(s_vec(igen,j),j=1,3)
157         endif
158c
159      elseif(itype.eq.0) then
160         if (oprint) write(*,31) 'PT'
161      endif
162c
16330    format(/,8x,a2,' fold Rotoinversion Operator at (',2(f9.6,','),f9.
164     &6,')')
16531    format(/,25x,a2,' fold Rotoinversion Operator')
166c
167c--> write out input generators
168c
169      if (oprint) call mprint(capr,3,4)
170c
171c--> get the current operator to do mults with, O(j)
172c
1735999  icnt2=1
1746000  jpos=(icnt2-1)*3+1
175      do 320 i=1,3
176         do 330 j=1,4
177            caps(i,j)=symops(jpos+(i-1),j)
178330      continue
179320   continue
180c
181c--> determine new rotational portion of operator [resop=capr*caps]
182c
183      do 130 i=1,3
184         do 140 j=1,3
185            sum=0.0d+00
186            do 150 k=1,3
187               sum = sum + capr(i,k)*caps(k,j)
188150         continue
189            resop(i,j)=sum
190140      continue
191130   continue
192c
193c--> compute new translational components [Ru+t]
194c
195      do 250 i=1,3
196         sum=0.0d+00
197         do 260 j=1,3
198            sum=sum+capr(i,j)*caps(j,4)
199260      continue
200         resop(i,4)=sum
201250   continue
202      do 270 i=1,3
203         resop(i,4)=resop(i,4)+capr(i,4)
204270   continue
205c
206c--> clean up translational components (put in standard convention)
207c
208      do 390 i=1,3
209         if(resop(i,4).ge.-tol.and.resop(i,4).le.tol) then
210            resop(i,4)=0.0d+00
211         elseif(resop(i,4).ge.1.0d+00-tol.and.resop(i,4).le.
212     &          1.0d+00+tol) then
213            resop(i,4)=1.0d+00
214         endif
215390   continue
216c
217c
218      do 400 i=1,3
219         if(resop(i,4).lt.0.0d+00) then
220            resop(i,4)=resop(i,4)+1.00000000
221         elseif(resop(i,4).ge.1.0d+00) then
222            resop(i,4)=resop(i,4)-1.00000000
223         endif
224400   continue
225c
226c---> Set translation components 1/3 and 2/3 of the generator
227c      products to double precision.
228c      Otherwise get 1/3=0.33333334 and 2/3=0.66666666.
229c      (A similar loop for the generators appears in "spgen.f")
230c
231      do 247 i=1,3
232         if(resop(i,4).le..6667.and.resop(i,4).ge..6666) then
233             resop(i,4)=2.0d+00/3.0d+00
234         elseif(resop(i,4).le..3334.and.resop(i,4).ge..3333) then
235             resop(i,4)=1.0d+00/3.0d+00
236         endif
237247   continue
238c
239c--> determine what the newly generated operator is
240c
241c**********************************************************************
242c     Move this after screens? (rotoop(j) is not changed here!)
243c**********************************************************************
244      do 160 i=1,3
245         do 170 j=1,3
246            detres(i,j)=resop(i,j)
247170      continue
248160   continue
249c
250c--> compute trace and determinant of the new operator
251c
252      trace=0.0d+00
253      do 180 i=1,3
254         trace=trace+resop(i,i)
255180   continue
256*
257      det = deter3(detres)
258**      call ludcmp(detres,3,3,indx,det)
259**      do 190 i=1,3
260**         det=det*detres(i,i)
261**190   continue
262c
263c-->  Check whether {RS|Ru+t} is identical to any of the previous
264c       operators stored in the "symops" matrix.
265c     The computation of Gi^2 (i>1) produces identical operators
266c       in groups #75-#230 (e.g. group #75(P4)).  Exclusion of Gi^2
267c       (i>1) from the algorithm fixes the problem for groups #75-#194,
268c       but it results in the loss of many operators from the cubic
269c       groups (e.g. group #195(P23)).
270c
271      irow=0
272      ichecks=1
273191   isame=0
274      do 192 i=1,3
275          do 193 j=1,4
276              if(resop(i,j).ge.symops(irow+i,j)-tol.and.resop(i,j).le.sy
277     &mops(irow+i,j)+tol) then
278                isame=isame+1
279              endif
280193       continue
281192    continue
282      if(isame.eq.12.and.jpos.lt.ipos) then
283          icnt2=icnt2+1
284          goto 6000
285      elseif(isame.eq.12.and.jpos.eq.ipos.and.igen.lt.ngen) then
286          icnt1=nops+1
287          igen=igen+1
288          goto 5000
289      elseif(isame.eq.12.and.jpos.eq.ipos.and.igen.eq.ngen) then
290          goto 223
291      elseif(ichecks.lt.nops) then
292          irow=irow+3
293          ichecks=ichecks+1
294          goto 191
295      endif
296c
297c--> add the new operator to the bottom of the symop list if it is not
298c    the identity element, increment pointer containing the number of
299c    operations currently in the symops list
300c
301c**********************************************************************
302c  Add the next screen since "idint" merely extracts the integer part
303c  of a real number. (E.g. idint(2.999999998912563)=2).
304c  Note: The addition is unnecessary for all 230 space groups and
305c        point groups 1-13. Point group #14 (D5d), where C5^5=E is
306c        computed, is the first case where "idint" causes problems.
307c                           DGC 2-25-94
308c**********************************************************************
309      if(trace.ge.3.00d+00-tol.and.trace.le.3.00d+00+tol) then
310        trace=3.00d+00
311      endif
312      itrace=idint(trace)
313      if(itrace.ne.3) then
314         isym=nops*3
315         do 350 i=1,3
316            do 360 j=1,4
317               symops(isym+i,j)=resop(i,j)
318360         continue
319350      continue
320         nops=nops+1
321c
322c--> assign labels to the rotational portion of the generated
323c    symmetry operation.
324c
325c**********************************************************************
326c  Screen out all point groups from the next calc   DGC 3/10/94
327c**********************************************************************
328c
329         if(det.lt.0.0d+00.and.itype.ne.0) then
330            rotoop(nops)=kneg(itrace)
331         elseif(itype.ne.0) then
332            rotoop(nops)=kpos(itrace)
333         endif
334      endif
335c
336c --> debug prints (print all matrices and new mats determinants)
337c
338c      write(*,12) 'trace = ', trace
339c      write(*,12) 'det   = ', det
340c      call mprint(resop,3,4)
341c----------------------------------------------------------------------
342c--> get ready to do next pass
343c
344c    If block:
345c       If ipos.lt.jpos then there are more operators in the list symops
346c      that need to be multiplied by the current generator.
347c
348c--->  If an identity was generated before exhausting all possible
349c     multipliers in symops, then get a new multiplier from symops.
350c      If an identity was generated after exhausting all possible
351c     multipliers in symops, then get a new generator if there is
352c     still one present in the generator list, and quit the algorithm
353c     if there aren't any generators left.
354c----------------------------------------------------------------------
355      if(itrace.eq.3.and.jpos.lt.ipos) then
356         icnt2=icnt2+1
357         goto 6000
358      elseif(itrace.eq.3.and.jpos.eq.ipos.and.igen.lt.ngen) then
359         icnt1=nops+1
360         igen=igen+1
361         goto 5000
362      elseif(itrace.eq.3.and.jpos.eq.ipos.and.igen.eq.ngen) then
363         goto 223
364      endif
365c
366c---------------------------------------------------------------------
367c      A new operator has been generated at this point, and if the
368c     list of possible multipliers in the list symops has not been
369c     exhausted, then get a new multiplier. If the list has been
370c     exhausted (jpos.eq.ipos), then get a new generator unless the
371c     space group has only one generator (ngen.eq.1).
372c      To compute terms such as G3^2G1, G3^2G1^1,...G3^2G2^2G1^2,
373c     the sequence jpos = 1,4,7,...,ipos of multipliers from symops
374c     is repeated - but this time they are left-multiplied by G3^2
375c     instead of G3. The indicator "isquare" = 0 or 1 when the left-
376c     multiplier is G3 or G3^2, rspt. At the end of the sequence
377c     having G1 as the left-multiplier, the matrix R (capr) is set
378c     equal to Gi^2.
379c---------------------------------------------------------------------
380      if(jpos.lt.ipos) then
381         icnt2=icnt2+1
382         goto 6000
383      elseif(igen.eq.1.and.igen.lt.ngen) then
384         igen=igen+1
385         icnt1=nops+1
386         goto 5000
387      elseif(ngen.eq.1) then
388         goto 223
389      elseif(igen.eq.ngen.and.isquare.eq.1) then
390         goto 223
391      elseif(isquare.eq.1) then
392         isquare=0
393         igen=igen+1
394         icnt1=nops+1
395         goto 5000
396      endif
397      do 201 i=1,3
398         do 221 j=1,4
399            capr(i,j)=symops(isym+i,j)
400221      continue
401201   continue
402      isquare=1
403      goto 5999
404c
40512    format(a8,f10.6)
406c
407c--> Determine if this is a centered space group, if so add the centering
408c    vectors to those generated for the point (0,0,0), ie to the current list
409c    of symops
410c
411c--> fill in roto-inversion operators (duplicate from top part of symops
412c    to the bottom) then add centering vectors to the 4th col. of symops
413c
414c
415c
416
417223   if(numvec.gt.0) then
418c------> fill in rotoinversion ops
419         icnt=nops*3
420         do 830 ij=1,numvec
421c------> Fill in the identity first
422            do 827 kiki=1,3
423               do 828 ikik=1,3
424                   symops(icnt+kiki,ikik)=0
425828            continue
426827         continue
427            do 829 kiki=1,3
428               symops(icnt+kiki,kiki)=1
429829         continue
430c---> label to rotop for identity
431            rotoop(nops*ij+ij)=' 1'
432            icnt=icnt+3
433            do 832 ktop=1,nops*3
434               icnt=icnt+1
435               do 834 kcol=1,3
436                  symops(icnt,kcol)=symops(ktop,kcol)
437834            continue
438832         continue
439830      continue
440c------> add centering vectors to col. 4
441         nops1=nops
442         do 849 i=1,numvec
443           isym=nops*3
444c------> Add cntvec to the identity first
445             ibot=isym+3*(i-1)
446             do 853 j=1,3
447                symops(ibot+j,4)=cntvec(i,j)
448853          continue
449            do 850 iold=1,nops1
450               itop=(iold-1)*3
451               ibot=isym+itop +3*i
452               do 860 j=1,3
453                  symops(ibot+j,4)=symops(itop+j,4)+cntvec(i,j)
454860            continue
455               rotoop(nops+i+iold)=rotoop(iold)
456850         continue
457            nops=nops1*(i+1)
458849      continue
459      endif
460c
461c--> add labels to centered symops
46225    format(a14,f10.6)
463      nops=nops+numvec
464c
465c--> print the matrix reps in operator form, with labels
466      if(oprint) then
467         write(*,1423) nops
468         call opprint(symops,rotoop,maxops,nops,itype)
469      endif
470c
471c dgc --> decenter if necessary
472      if(numvec.gt.0) then
473         !write(*,1420)
474         use_primitive = geom_use_primitive(geom)
475         if (use_primitive) then
476            write(*,1421) group_name
477            call dctr(symops,nops1,numgrp,group_name,numvec,cntvec,
478     >                numset,geom)
479            write(*,1424)
480            write(*,1425)
481            write(*,1426)
482            write(*,1427)
483            write(*,1428)
484            write(*,1429) nops+1,nops1+1
485            nops=nops1
486            if (oprint) then
487               write(*,1431)
488               call opprint(symops,rotoop,maxops,nops,itype)
489            end if
490            !write(*,1430)
491         else
492            write(*,1422) group_name
493         end if
494      endif
495c
4961420  format(/' Primitive cell exists. ')
4971421  format(/' Primitive cell exists. Converting the symmetry',
498     & ' operators and the unit cell to primitive cell.'/,
499     & ' To not do this conversion',
500     & ' add the "conventional" keyword to the symmetry input, e.g.'//,
501     & ' geometry'/," ..."/,' symmetry ',A,' conventional'/,
502     & ' ...'/,' end'/)
5031422  format(/' Primitve cell exists, but not converting to it.',
504     & ' Using the conventional cell',
505     & ' instead.'/,' To turn this conversion on',
506     & ' add the "primitive" keyword to the symmetry input, e.g.'//,
507     & ' geometry'/," ..."/,' symmetry ',A,' primitive'/,
508     & ' ...'/,' end'/)
5091423  format(//'The ',i3,' symmetry operators (excl. E)',
510     &' are as follows:'/)
5111424  format(/'After de-centering the following are redefined:'/)
5121425  format(10X,' (1) lattice parameters a,b,c,alpha,beta, and gamma')
5131426  format(10X,' (2) lattice vectors a1,a2,a3')
5141427  format(10X,' (3) reciprocal lattice vectors b1,b2,b3')
5151428  format(10X,' (4) fractional coordinates')
5161429  format(/i3,' operators converted to ',i3,' symmetry operators.')
5171430  format(/'DEBUG:primitive cell exists, but dctr was not called.'/)
5181431  format(//i3,'The symmetry operators (excl. E) in the de-centered',
519     &' basis are as follows:'/)
520c
521c rjh hack to fix C1
522      if (numgrp .eq. 1) nops = 0
523c
524      end
525
526