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