1 SUBROUTINE SYMPOP(H,I,ISKIP,DELDIP) 2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 3 INCLUDE 'SIZES' 4 DIMENSION H(*), DELDIP(3,*) 5 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 6 COMMON /SYMOPC/ ISYMT(6) 7 CHARACTER*10 ISYMT 8 COMMON /ATOMS/ COORD, NATOMS, NVAR 9 DO 10 J = 1, NSYM 10 IF (IPO(I,J).LT.I) THEN 11 CALL SYMH(H, DELDIP, I, J) 12 ISKIP=3 13C atom ipo(i,j) is suitable for transition dipole calc'n 14C 15 K=I*3-2 16C 17C INSERT DELDIP ROTATION HERE 18C 19 GOTO 20 20 ENDIF 21 10 CONTINUE 22 ISKIP=0 23 20 CONTINUE 24 RETURN 25 END 26 SUBROUTINE SYMR 27 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 28 INCLUDE 'SIZES' 29 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 30 COMMON /SYMOPC/ ISYMT(6) 31 CHARACTER*10 ISYMT 32***************************************************************** 33* 34* ON INPUT NONE 35* 36* ON OUTPUT R = SYMMETRY OPERATIONS 37* IPO = PERMUTATION OPERATOR FOR SYMMETRY OPERATIONS 38* NSYM = NUMBER OF SYMMETRY OPERATIONS 39* NENT = NUMBER OF SYMMETRY OPERATIONS ENTERED 40* 41***************************************************************** 42C 43C A SUBROUTINE THAT WILL READ IN THE PRIMATIVE SYMMETRY OPERATIONS 44C AND DETERMINE IF THEY ARE VALID FOR THIS MOLECULE. THIS 45C INFORMATION IS THEN EXPANDED TO THE COMPLETE SET AND USED FOR 46C SYMMATRIZING THE HESSIAN. 47C 48C THE CORRECT FORMAT FOR DESCRIBING A SYMMETRY FUNCTION IS: 49C LABEL IFNCN AXIS 50C 51C WHERE: 52C LABEL - MUST BE INCLUDED AND IS THE LABEL THAT WILL 53C BE USED TO IDENTIFY THAT FUNCTION 54C IFNCN - THE NUMBER OF THE SYMMETRY FUNCTION TO BE USED: 55C 0 - INVERSION OPERATOR 56C 1 - REFLECTION PLANE PERPENDICULAR TO THE AXIS 57C 2-X - A C(N) AXIS 58C -3-X - A S(N) AXIS 59C AXIS - THE AXIS FOR THE OPERATION. MAY BE SPECIFIED AS: 60C X, Y, Z COORDINATES -> MUST USE 3 VALUES. AT LEAST ONE 61C MUST BE NON-INTEGER OR TWO MUST BE IDENTICAL 62C AN ATOM LIST -> THE COORDINATES OF THE ATOMS LISTED WILL 63C BE SUMMED TO GENERATE THE AXIS. 64C A -B -> THE VECTOR FROM ATOM B TO ATOM A WILL BE USED AS 65C THE AXIS. 66C 67C A MAXIMUM OF 6 SYMMETRY OPERATIONS CAN BE INPUTTED. THESE SHOULD BE THE 68C UNIQUE GENERATING FUNCTIONS FROM WHICH ALL THE OPERATIONS OF THE GROUP 69C CAN BE CONSTRUCTED. E.G. ONLY C5 NEEDS TO BE SPECIFIED SINCE C5(2) 70C THROUGH C5(4) CAN BE GENERATED FROM THIS SINGLE FUNCTION. 71C 72C A MAXIMUM OF 8 UNIQUE ATOMS CAN BE USED TO SPECIFY AN AXIS. 73C 74C THE E FUNCTION IS BY DEFAULT THE FIRST SYMMETRY FUNCTION. THIS FUNCTION 75C NEVER NEEDS TO BE EXPLICTLY INCLUDED IN YOUR LIST. IT CANNOT BE 76C ENTERED. 77C 78C IF YOU ENTER A GIVEN SYMMETRY FUNCTION MORE THAN ONCE, ONLY THE FIRST 79C OCCURANCE WILL BE USED. ALL DUPLICATES WILL BE DELETED. 80C 81 DIMENSION TEMP(9), TEMP2(9), ISTART(7) 82 INTEGER ITEMP(9) 83 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 84 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 85 2 NCLOSE,NOPEN,NDUMY,FRACT 86 COMMON /COORD / COORD(3,NUMATM) 87 COMMON /KEYWRD/ KEYWRD 88 CHARACTER KEYWRD*241 89 CHARACTER*80 LINE 90 LOGICAL LEADSP, PROB, ALLINT 91C THE NEXT PARAMETERS ARE THE MAX NUMBER OF SYMM FUNCTIONS, THE 92C MAX NUMBER OF SYMM FUNCTIONS TO READ IN, AND THE 93C TOLERENCE USED TO DETERMINE IF TWO FUNCTIONS ARE IDENTICAL 94 PARAMETER (MAXFUN=120) 95 PARAMETER (MAXENT=6) 96 PARAMETER (TOL =1D-3) 97C 98C 99C Variables used: (n represents the number of atomic centers) 100C TEMP(9), TEMP2(9): Temporary matricies used to hold small parts 101C larger matricies for specific matrix operations. 102C 103C For the next two items, the last index represents the symmetry 104C operation number. 105C R(9,*): The 9 elements of each record are a packed 3 by 3 106C array of a given symmetry operations. 107C IPO(n,*): A vector that contains the symmetry mapping of atomic ce 108C 109 PROB = .FALSE. 110C Get the symmetry functions: (NOTE: THE FIRST IS ALWAY E) 111 R(1,1) = 1.D0 112 R(2,1) = 0.D0 113 R(3,1) = 0.D0 114 R(4,1) = 0.D0 115 R(5,1) = 1.D0 116 R(6,1) = 0.D0 117 R(7,1) = 0.D0 118 R(8,1) = 0.D0 119 R(9,1) = 1.D0 120C 121 X = 0.D0 122 Y = 0.D0 123 Z = 0.D0 124 DO 10 I=1,NUMAT 125 X = X + COORD(1,I) 126 Y = Y + COORD(2,I) 127 Z = Z + COORD(3,I) 128 IPO(I,1) = I 129 10 CONTINUE 130 XA = X/FLOAT(NUMAT) 131 YA = Y/FLOAT(NUMAT) 132 ZA = Z/FLOAT(NUMAT) 133 DO 20 I=1,NUMAT 134 COORD(1,I) = -XA + COORD(1,I) 135 COORD(2,I) = -YA + COORD(2,I) 136 COORD(3,I) = -ZA + COORD(3,I) 137 20 CONTINUE 138 WRITE(6,'(/'' SYMMETRY OPERATIONS USED FOR SYMMETRIZING'', 139 1'' THE HESSIAN'')') 140 WRITE(6,'(/,'' OPERATOR TYPE AXIS DEFINITION '')') 141C 142 NENT = 1 143 NSYM = 0 144 ISYMT(1)='E' 145 30 NSYM = NSYM + 1 146 READ(5,'(A)',END=120,ERR=120)LINE 147 LEADSP=.TRUE. 148 NVALUE=0 149 ALLINT=.TRUE. 150 DO 40 I=1,80 151 IF (LEADSP.AND.LINE(I:I).NE.' ') THEN 152 NVALUE=NVALUE+1 153 ISTART(NVALUE)=I 154 ENDIF 155 LEADSP=(LINE(I:I).EQ.' ') 156 40 CONTINUE 157 IF (NVALUE.EQ.0) GOTO 120 158 IF (NVALUE.EQ.1) THEN 159 WRITE(6,200) 160 200 FORMAT(' NOT A VALID LINE. ONLY HAS ONE ENTRY') 161 PROB = .TRUE. 162 GOTO 120 163 ENDIF 164 ISYMT(1+NENT)=LINE(ISTART(1):ISTART(2)-1) 165 DO 50 I=2,NVALUE 166 TEMP(I-1)=READA(LINE,ISTART(I)) 167 ITEMP(I-1)=NINT(TEMP(I-1)) 168 IF((ABS(ITEMP(I-1)-TEMP(I-1)).GT.TOL).AND.(I.NE.2)) 169 + ALLINT=.FALSE. 170 50 CONTINUE 171 IF (ALLINT) THEN 172 WRITE(6,210)ISYMT(1+NENT),(ITEMP(I),I=1,NVALUE-1) 173 210 FORMAT(X,A10,I7,8I7) 174 ELSE 175 WRITE(6,220)ISYMT(1+NENT),ITEMP(1),(TEMP(I),I=2,NVALUE-1) 176 220 FORMAT(X,A10,I7,8F7.3) 177 ENDIF 178 SIGMA = 1 179 IF (ITEMP(1) .LE. -3) SIGMA = -1 180 TEMP(1) = ABS( TEMP(1)) 181 ITEMP(1)= ABS(ITEMP(1)) 182 IF (ABS(ITEMP(1)-TEMP(1)) .GE. TOL) THEN 183 WRITE(6,230) 184 230 FORMAT(' THE SYMMETRY FUNCTION MUST BE INTEGER') 185 PROB = .TRUE. 186 GOTO 120 187 ENDIF 188 IF (ITEMP(1) .EQ. 0) THEN 189C WITH INVERSION, THE AXIS IS UNIMPORTANT 190 R(1,1+NENT) = -1.D0 191 R(2,1+NENT) = 0.D0 192 R(3,1+NENT) = 0.D0 193 R(4,1+NENT) = 0.D0 194 R(5,1+NENT) = -1.D0 195 R(6,1+NENT) = 0.D0 196 R(7,1+NENT) = 0.D0 197 R(8,1+NENT) = 0.D0 198 R(9,1+NENT) = -1.D0 199 GOTO 70 200 ENDIF 201C WITH ANYTHING ELSE, THE AXIS MUST BE DETERMINED. IF NO AXIS IS DEFINED 202C FLAG IT AS A PROBLEM 203 IF (NVALUE .EQ. 2) THEN 204 PROB = .TRUE. 205 WRITE(6,240) NSYM 206 240 FORMAT(' NO AXIS INFORMATION WAS ENTERED FOR FUNCTION',I2) 207 GOTO 120 208 ENDIF 209 IF ((NVALUE .EQ. 5).AND.((TEMP(2).EQ.TEMP(3)) 210 + .OR.(TEMP(2).EQ.TEMP(4)).OR. 211 + (TEMP(3).EQ.TEMP(4)).OR.(.NOT. ALLINT).OR. 212 + (ABS(ITEMP(2)).LT. 1).OR.(ABS(ITEMP(2)).GT.NUMAT) .OR. 213 + (ABS(ITEMP(3)).LT. 1).OR.(ABS(ITEMP(3)).GT.NUMAT) .OR. 214 + (ABS(ITEMP(4)).LT. 1).OR.(ABS(ITEMP(4)).GT.NUMAT))) THEN 215C IT APPEARS TO BE XYZ INPUT 216 X = TEMP(2) 217 Y = TEMP(3) 218 Z = TEMP(4) 219 ELSE 220C APPEARS TO BE ATOM NUMBER INPUT 221 IF (.NOT. ALLINT) THEN 222 PROB = .TRUE. 223 WRITE(6,250) 224 250 FORMAT(' YOU MUST HAVE ALL INTEGER INPUT WHEN NOT', 225 + ' USING XYZ INPUT') 226 GOTO 120 227 ENDIF 228 X = 0.D0 229 Y = 0.D0 230 Z = 0.D0 231 DO 60 I = 2, NVALUE-1 232 IF ((ABS(ITEMP(I)).LT.1).OR.(ABS(ITEMP(I)).GT.NUMAT)) THEN 233 WRITE(6,260)ITEMP(I) 234 260 FORMAT(' ATOM NUMBER',I3,' IS OUT OF RANGE') 235 PROB=.TRUE. 236 ENDIF 237 X = X + ITEMP(I)/ABS(ITEMP(I))*COORD(1,ABS(ITEMP(I))) 238 Y = Y + ITEMP(I)/ABS(ITEMP(I))*COORD(2,ABS(ITEMP(I))) 239 Z = Z + ITEMP(I)/ABS(ITEMP(I))*COORD(3,ABS(ITEMP(I))) 240 60 CONTINUE 241 ENDIF 242C 243C TIME TO DECIPHER THE SYMMETRY FUNCTION 244C 245 IF (ITEMP(1) .GT. 10) THEN 246 WRITE(6,270) 247 270 FORMAT(' A C-10 AXIS IS THE HIGHEST THAT CAN BE SPECIFIED') 248 PROB = .TRUE. 249 GOTO 120 250 ENDIF 251 ROT=4.D0*ASIN(1.D0)/ITEMP(1) 252 IF(ITEMP(1).EQ. 1) SIGMA = -1 253C 254C First, construct the matrix defining the rotation axis 255 XY=X**2+Y**2 256 RA=SQRT(XY+Z**2) 257 IF (RA.LT. TOL) THEN 258 PROB = .TRUE. 259 WRITE(6,280) 260 280 FORMAT(' YOUR VECTOR AXIS MUST HAVE A NON-ZERO LENGTH ') 261 GOTO 120 262 ENDIF 263 XY=SQRT(XY) 264 IF (XY.GT.1.D-10) THEN 265 CA=Y/XY 266 CB=Z/RA 267 SA=X/XY 268 SB=XY/RA 269 ELSEIF (Z.GT. 0.D0) THEN 270 CA=1.D0 271 CB=1.D0 272 SA=0.D0 273 SB=0.D0 274 ELSE 275 CA=-1.D0 276 CB=-1.D0 277 SA=0.D0 278 SB=0.D0 279 ENDIF 280C GENERATE THE MATRIX ELEMENTS BY DOING THE EULER TRANSFORM 281 TEMP( 1)=CA 282 TEMP( 2)=-SA 283 TEMP( 3)=0.D0 284 TEMP( 4)=SA*CB 285 TEMP( 5)=CA*CB 286 TEMP( 6)=-SB 287 TEMP( 7)=SA*SB 288 TEMP( 8)=CA*SB 289 TEMP( 9)=CB 290C 291C 292 CA = DCOS(ROT) 293 SA = DSIN(ROT) 294C 295C The construct the actual R matrix to be used 296C 297C 298 TEMP2(1) = CA 299 TEMP2(2) = SA 300 TEMP2(3) = 0.D0 301 TEMP2(4) = -SA 302 TEMP2(5) = CA 303 TEMP2(6) = 0.D0 304 TEMP2(7) = 0.D0 305 TEMP2(8) = 0.D0 306 TEMP2(9) = SIGMA 307 CALL MAT33(TEMP, TEMP2, R(1,1+NENT)) 308C 309C Now, verify that this is a unique and valid function 310 70 CONTINUE 311C 312 RES = 10.D0 313 DO 90 I = 2, NENT 314 RESO = 0.D0 315 DO 80 J = 1, 9 316 80 RESO= ABS( R(J,I) - R(J,1+NENT)) + RESO 317 RES = MIN(RES, RESO) 318 90 CONTINUE 319 IF (RES .LT. TOL) THEN 320C THIS IS NOT VALID FUNCTION 321 WRITE(6,290) 322 290 FORMAT(' THIS FUNCTION IS IDENTICAL TO AN EARLIER ONE') 323 GOTO 120 324 ENDIF 325C NOW, TO CALCULATE THE IPO OF THIS FUNCTION 326 NENT = 1+NENT 327 N = NENT 328C Now, to initialize IPO(n) and 329C Perform R on each atomic center and determine where it maps to. 330 DO 110 I = 1, NUMAT 331 X=COORD(1,I)*R(1,N) + COORD(2,I)*R(2,N) + COORD(3,I)*R(3,N) 332 Y=COORD(1,I)*R(4,N) + COORD(2,I)*R(5,N) + COORD(3,I)*R(6,N) 333 Z=COORD(1,I)*R(7,N) + COORD(2,I)*R(8,N) + COORD(3,I)*R(9,N) 334 IPO(I,N) = 0 335 DO 100 J = 1, NUMAT 336 DIST=ABS(X-COORD(1,J))+ABS(Y-COORD(2,J))+ABS(Z-COORD(3,J)) 337 IF (DIST .LT. 5.D-2) THEN 338 IF (IPO(I,N) .EQ. 0) THEN 339 IPO(I,N) = J 340 ELSE 341 WRITE(6,300) 342 PROB = .TRUE. 343 GOTO 120 344 300 FORMAT(' ONE ATOM MAPS ONTO TWO DIFFERENT ATOMIC C', 345 1'ENTERS') 346 ENDIF 347 ENDIF 348 100 CONTINUE 349 IF (IPO(I,N) .EQ. 0) THEN 350 WRITE(6,310) 351 310 FORMAT(' ONE ATOM MAPS ONTO NO OTHER ATOM ') 352 PROB = .TRUE. 353 GOTO 120 354 ENDIF 355 110 CONTINUE 356C 357C IF THIS POINT IS REACHED, THE FUNCTION IS VALID 358C CHECK IF THE R MATRIX SHOULD BE PRINTED 359C 360 IF (INDEX(KEYWRD,' RMAT') .NE. 0) THEN 361 WRITE(6,320)(R(I,N),I=1,3) 362 WRITE(6,330)N,(R(I,N),I=4,6) 363 WRITE(6,340)(R(I,N),I=7,9) 364 320 FORMAT(/,10X,'| ',3F10.6,' |') 365 330 FORMAT(I5,' = | ',3F10.6,' |') 366 340 FORMAT(10X,'| ',3F10.6,' |',/) 367 ENDIF 368C 369 120 IF((NVALUE.NE.0).AND.(NSYM.LT.MAXENT)) GOTO 30 370C 371C If a problem exists. Stop the program. 372C 373 IF (PROB) THEN 374 CLOSE (6) 375 STOP 'PROBLEM IN SYMR' 376 ENDIF 377C 378C NOW, ALL USER FUNCTIONS ARE IN WITH NO ERRORS (JUST ELIMINATION OF DUPS) 379C 380 IF(NVALUE.NE.0) READ(5,'(A)',END=130)LINE 381 130 CONTINUE 382 NSYM = NENT 383C 384C NEXT, EXPAND THE EXISTING OPERATORS TO THE FULL SET 385C 386 CALL SYMP 387C 388 RETURN 389 END 390 SUBROUTINE SYMH(H, DIP, I, N) 391 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 392 INCLUDE 'SIZES' 393 DIMENSION H(*), DIP(3,*) 394 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 395 COMMON /SYMOPC/ ISYMT(6) 396 CHARACTER*10 ISYMT 397***************************************************************** 398* 399C INPUT: H() A packed lower triangular hessian 400C DIP(,) A MATRIX OF DIPOLE TENSORS TO BE SYMM 401C R(,) A matrix of symmetry operations 402C IPO(,) A matrix of atomic mapping according to R 403C I The atom (row and column) to add to H() 404C N The symmetry operation to use to generate I 405C 406C OUTPUT: H() A packed lower triangular Hessian with information 407C about atom I added 408C DIP(,) A MATRIX OF DIPOLE TENSORS THAT HAVE BEEN SYMM 409C 410***************************************************************** 411C 412C 413C This subroutine will add all necessary information to the Hessian con 414C atom I. Since the Hessian is a packed lower half triangle, the exi 415C information for atom pair (K,L) where K,L < I is fully known, (K > 416C L < I) or (vice versa) is half known, K,L > I is completely unknown 417C Therefore, start in unknown region and make it half known. Double 418C known values, and move in the diagonal element at full strength. 419C 420C 421 DIMENSION TEMP(9), TEMP2(9) 422 COMMON /FOKMAT/ HA(MPACK*2) 423 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 424 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 425 2 NCLOSE,NOPEN,NDUMY,FRACT 426 COMMON /COORD / COORD(3,NUMATM) 427C 428C 429C 430C Variables used: (n represents the number of atomic centers) 431C H(3n,3n): Input/output matrix. It is a packed lower half triangu 432C matrix. Commonly, the Hessian. 433C TEMP(9), TEMP2(9): Temporary matricies used to hold small parts 434C larger matricies for specific matrix operations. 435C 436C For the next two items, the last indicy represents the symmetry 437C operation number. 438C R(14,*): The first 9 elements of each record are a packed 3 by 3 439C array of a given symmetry operations. Elements 10 - 14 are t 440C users input describing the symmetry operation. 441C IPO(n,*): A vector that contains the symmetry mapping of atomic c 442C 443C 444 K = IPO(I,N) 445C 446C Now, to climb up the matrix 447 DO 10 J = NUMAT, I+1, -1 448 L = IPO(J,N) 449C 450C Now, to actually perform R H R 451C 452C Do this multiplication in a 3 by 3 block at a time. Store H(i,j) in 453C H( IPO(I,N), IPO(J,N)) 454C 455 IF (K .GT. L) THEN 456 IEL33 = (3*K*(3*K-1))/2 + 3*L 457 TEMP(9) = 0.5D0 * H(IEL33) 458 TEMP(8) = 0.5D0 * H(IEL33-1) 459 TEMP(7) = 0.5D0 * H(IEL33-2) 460 TEMP(6) = 0.5D0 * H(IEL33-K*3+1) 461 TEMP(5) = 0.5D0 * H(IEL33-K*3) 462 TEMP(4) = 0.5D0 * H(IEL33-K*3-1) 463 TEMP(3) = 0.5D0 * H(IEL33-6*K+3) 464 TEMP(2) = 0.5D0 * H(IEL33-6*K+2) 465 TEMP(1) = 0.5D0 * H(IEL33-6*K+1) 466 ELSE 467 IEL33 = (3*L*(3*L-1))/2 + 3*K 468 FACT = 1.0D0 469 IF (L .LT. I) FACT = 0.5D0 470 TEMP(9) = FACT * H(IEL33) 471 TEMP(6) = FACT * H(IEL33-1) 472 TEMP(3) = FACT * H(IEL33-2) 473 TEMP(8) = FACT * H(IEL33-L*3+1) 474 TEMP(5) = FACT * H(IEL33-L*3) 475 TEMP(2) = FACT * H(IEL33-L*3-1) 476 TEMP(7) = FACT * H(IEL33-6*L+3) 477 TEMP(4) = FACT * H(IEL33-6*L+2) 478 TEMP(1) = FACT * H(IEL33-6*L+1) 479 ENDIF 480C 481 CALL MAT33(R(1,N), TEMP, TEMP2) 482C 483 IEL33 = J*3*(J*3-1)/2 + I*3 484 H(IEL33) = TEMP2(9) 485 H(IEL33-3*J+1) = TEMP2(8) 486 H(IEL33-6*J+3) = TEMP2(7) 487 H(IEL33-1) = TEMP2(6) 488 H(IEL33-3*J) = TEMP2(5) 489 H(IEL33-6*J+2) = TEMP2(4) 490 H(IEL33-2) = TEMP2(3) 491 H(IEL33-3*J-1) = TEMP2(2) 492 H(IEL33-6*J+1) = TEMP2(1) 493 10 CONTINUE 494C 495C Now, to do the diagonal term 496C 497 IEL33 = (3*K*(3*K+1))/2 498 TEMP(9) = 0.5D0 * H(IEL33) 499 TEMP(8) = 0.5D0 * H(IEL33-1) 500 TEMP(7) = 0.5D0 * H(IEL33-2) 501 TEMP(6) = TEMP(8) 502 TEMP(5) = 0.5D0 * H(IEL33-K*3) 503 TEMP(4) = 0.5D0 * H(IEL33-K*3-1) 504 TEMP(3) = TEMP(7) 505 TEMP(2) = TEMP(4) 506 TEMP(1) = 0.5D0 * H(IEL33-6*K+1) 507C 508 CALL MAT33(R(1,N), TEMP, TEMP2) 509C 510 IEL33 = I*3*(I*3+1)/2 511 H(IEL33) = TEMP2(9) 512 H(IEL33-1) = TEMP2(8) 513 H(IEL33-2) = TEMP2(7) 514 H(IEL33-I*3) = TEMP2(5) 515 H(IEL33-I*3-1) = TEMP2(4) 516 H(IEL33-6*I+1) = TEMP2(1) 517C 518C NOW, TO ROTATE THE DIPOLE TENSOR TERM 519C 520 TEMP(9) = DIP(3,K*3 ) 521 TEMP(8) = DIP(2,K*3 ) 522 TEMP(7) = DIP(1,K*3 ) 523 TEMP(6) = DIP(3,K*3-1) 524 TEMP(5) = DIP(2,K*3-1) 525 TEMP(4) = DIP(1,K*3-1) 526 TEMP(3) = DIP(3,K*3-2) 527 TEMP(2) = DIP(2,K*3-2) 528 TEMP(1) = DIP(1,K*3-2) 529C 530 CALL MAT33(R(1,N), TEMP, TEMP2) 531C 532 DIP(3,I*3 ) = TEMP2(9) 533 DIP(2,I*3 ) = TEMP2(8) 534 DIP(1,I*3 ) = TEMP2(7) 535 DIP(3,I*3-1) = TEMP2(6) 536 DIP(2,I*3-1) = TEMP2(5) 537 DIP(1,I*3-1) = TEMP2(4) 538 DIP(3,I*3-2) = TEMP2(3) 539 DIP(2,I*3-2) = TEMP2(2) 540 DIP(1,I*3-2) = TEMP2(1) 541C 542C 543C Now, to double all existing values going across 544 ISTART = (I-1)*3*((I-1)*3+1)/2+1 545 DO 20 J = ISTART, IEL33 546 H(J) = H(J) + H(J) 547 20 CONTINUE 548C Everything is now done for this symmetry element. 549C 550 RETURN 551 END 552 SUBROUTINE SYMA(E, V) 553 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 554 INCLUDE 'SIZES' 555 DIMENSION E(NUMAT*3), V(NUMAT*NUMAT*9) 556 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 557 COMMON /SYMOPC/ ISYMT(6) 558 CHARACTER*10 ISYMT 559********************************************************************* 560* 561* ON INPUT E = FREQUENCIES IN CM(-1) 562* V = EIGENVECTORS OF NORMAL MODES, NORMALIZED 563* R = SYMMETRY OPERATIONS 564* IPO = MAP OF ATOMS BEING MOVED 565* NSYM = NUMBER OF SYMMETRY OPERATION 566* 567********************************************************************* 568C 569C THIS SUBROUTINE DETERMINES THE SYMMETRY FUNCTION VALUE OF EACH 570C VIBRATIONAL MODE. IT DOES IT BY DOING <EV R EV> 571C 572 COMMON /COORD / COORD(3,NUMATM) 573 COMMON /KEYWRD/ KEYWRD 574 DIMENSION T1(MAXPAR), T2(MAXPAR,7) 575 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 576 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 577 2 NCLOSE,NOPEN,NDUMY,FRACT 578 CHARACTER KEYWRD*241 579C 580 TOL=1.D-3 581C 582 NVAR=NUMAT*3 583C 584C T1(NVAR) AND T2(NVAR,NSYM) ARE THE ONLY ADDITIONAL ARRAYS NEEDED. TH 585C ARE TEMPORARY ARRAYS. 586 DO 30 K = 0, NVAR-1 587 DO 20 N = 1, NENT 588 DO 10 I = 1, NUMAT 589C 590 J=IPO(I,N) 591 T1(I*3-2)=V(J*3-2+K*NVAR)*R(1,N)+ 592 1 V(J*3-1+K*NVAR)*R(4,N)+ 593 2 V(J*3 +K*NVAR)*R(7,N) 594 T1(I*3-1)=V(J*3-2+K*NVAR)*R(2,N)+ 595 1 V(J*3-1+K*NVAR)*R(5,N)+ 596 2 V(J*3 +K*NVAR)*R(8,N) 597 T1(I*3 )=V(J*3-2+K*NVAR)*R(3,N)+ 598 1 V(J*3-1+K*NVAR)*R(6,N)+ 599 2 V(J*3 +K*NVAR)*R(9,N) 600 10 CONTINUE 601 T2(K+1,N) = 0.0D0 602 DO 20 I = 1, NVAR 603 T2(K+1,N) = T2(K+1,N) + T1(I)*V(I+K*NVAR) 604 20 CONTINUE 605 30 CONTINUE 606 WRITE(6,100) 607 WRITE(6,'('' '',7A9)')(ISYMT(I),I=1,NENT) 608 100 FORMAT(' FREQ.',/,' NO. FREQ. CHARACTER TABLE ') 609 I=1 610 J=I+1 611 IF (INDEX(KEYWRD,' NODEGEN') .NE. 0) TOL = -1.D0 612 EREF = E(1) 613 110 IF(ABS((E(J)-EREF)) .LE. TOL) THEN 614 DO 120 K = 1, NENT 615 120 T2(I,K) = T2(I,K) + T2(J,K) 616 E(I) = (E(I) + E(J)) 617 J = J+1 618 ELSE 619 E(I)=E(I)/FLOAT(J-I) 620 WRITE(6,130)I,E(I),(T2(I,K),K=1,NENT) 621 I=J 622 J=J+1 623 EREF=E(I) 624 ENDIF 625 IF (J .LE. NVAR) GOTO 110 626 E(I)=E(I)/FLOAT(J-I) 627 WRITE(6,130)I,E(I),(T2(I,K),K=1,NENT) 628 130 FORMAT(I4,F9.3,3X,7F9.4) 629 END 630 SUBROUTINE SYMT(H, DELDIP) 631 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 632 INCLUDE 'SIZES' 633 DIMENSION H(*), DELDIP(3,*) 634 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 635 COMMON /SYMOPC/ ISYMT(6) 636 CHARACTER*10 ISYMT 637***************************************************************** 638* 639* ON INPUT H = HESSIAN MATRIX, PACKED LOWER HALF TRIANGLE 640* R = SYMMETRY OPERATIONS 641* IPO = MAP OF ATOMS MOVED 642* NSYM = NUMBER OF SYMMETRY OPERATIONS 643* 644* ON OUTPUT H = SYMMETRIZED HESSIAN MATRIX 645* 646***************************************************************** 647C A subroutine that will symmatrize the Hamiltonian, or other matrix 648C by successive application of group operations. The method used 649C is R H R added to HA then divided by the total number of symmetry 650C operations used. This in effects averages all the values in a 651C symmetry correct fashion. 652C 653 DIMENSION TEMP(9), TEMP2(9), DELTMP(3,MAXPAR) 654 COMMON /FOKMAT/ HA(MPACK*2) 655 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 656 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 657 2 NCLOSE,NOPEN,NDUMY,FRACT 658 COMMON /COORD / COORD(3,NUMATM) 659C 660C 661C Variables used: (n represents the number of atomic centers) 662C H(3n,3n): Input/output matrix. It is a packed lower half triangu 663C matrix. Commonly, the Hessian. 664C HA(3n,3n): An internal matrix used to sum the symatrized Hessian 665C NSYM: Input, the value of this symmetry operation. 666C TEMP(9), TEMP2(9): Temporary matricies used to hold small parts 667C larger matricies for specific matrix operations. 668C 669C For the next two items, the last indicy represents the symmetry 670C operation number. 671C IPO(n,*): A vector that contains the symmetry mapping of atomic c 672C 673C Skip this subroutine if NSYMM <= 0. This implies that only E is pre 674 IF (NSYM .LT. 2) RETURN 675C 676 DO 10 I=1,(3*NUMAT*(NUMAT*3+1))/2 677 10 HA(I)=0.D0 678C 679 DO 15 I = 1, NUMAT*3 680 DELTMP(1,I) = 0.D0 681 DELTMP(2,I) = 0.D0 682 15 DELTMP(3,I) = 0.D0 683C 684 DO 40 N = 1, NSYM 685C 686C Now, to actually perform R H R 687 DO 30 I = 1, NUMAT 688 DO 20 J = 1, I-1 689C 690C Do this multiplication in a 3 by 3 block at a time. Store H(i,j) in 691C HA( IPO(I,N), IPO(J,N)) or HS( IPO(I,N), IPO(J,N)) 692C 693 K = IPO(I,N) 694 L = IPO(J,N) 695 IF (K .GT. L) THEN 696 IEL33 = (3*K*(3*K-1))/2 + 3*L 697 TEMP(9) = H(IEL33) 698 TEMP(8) = H(IEL33-1) 699 TEMP(7) = H(IEL33-2) 700 TEMP(6) = H(IEL33-K*3+1) 701 TEMP(5) = H(IEL33-K*3) 702 TEMP(4) = H(IEL33-K*3-1) 703 TEMP(3) = H(IEL33-6*K+3) 704 TEMP(2) = H(IEL33-6*K+2) 705 TEMP(1) = H(IEL33-6*K+1) 706 ELSE 707 IEL33 = (3*L*(3*L-1))/2 + 3*K 708 TEMP(9) = H(IEL33) 709 TEMP(6) = H(IEL33-1) 710 TEMP(3) = H(IEL33-2) 711 TEMP(8) = H(IEL33-L*3+1) 712 TEMP(5) = H(IEL33-L*3) 713 TEMP(2) = H(IEL33-L*3-1) 714 TEMP(7) = H(IEL33-6*L+3) 715 TEMP(4) = H(IEL33-6*L+2) 716 TEMP(1) = H(IEL33-6*L+1) 717 ENDIF 718C 719 CALL MAT33(R(1,N), TEMP, TEMP2) 720C 721 IEL33 = I*3*(I*3-1)/2 + J*3 722 HA(IEL33) = TEMP2(9) + HA(IEL33) 723 HA(IEL33-1) = TEMP2(8) + HA(IEL33-1) 724 HA(IEL33-2) = TEMP2(7) + HA(IEL33-2) 725 HA(IEL33-I*3+1) = TEMP2(6) + HA(IEL33-I*3+1) 726 HA(IEL33-I*3) = TEMP2(5) + HA(IEL33-I*3) 727 HA(IEL33-I*3-1) = TEMP2(4) + HA(IEL33-I*3-1) 728 HA(IEL33-6*I+3) = TEMP2(3) + HA(IEL33-6*I+3) 729 HA(IEL33-6*I+2) = TEMP2(2) + HA(IEL33-6*I+2) 730 HA(IEL33-6*I+1) = TEMP2(1) + HA(IEL33-6*I+1) 731 20 CONTINUE 732 K = IPO(I,N) 733 IEL33 = (3*K*(3*K+1))/2 734 TEMP(9) = H(IEL33) 735 TEMP(8) = H(IEL33-1) 736 TEMP(7) = H(IEL33-2) 737 TEMP(6) = TEMP(8) 738 TEMP(5) = H(IEL33-K*3) 739 TEMP(4) = H(IEL33-K*3-1) 740 TEMP(3) = TEMP(7) 741 TEMP(2) = TEMP(4) 742 TEMP(1) = H(IEL33-6*K+1) 743C 744 CALL MAT33(R(1,N), TEMP, TEMP2) 745C 746 IEL33 = I*3*(I*3+1)/2 747 HA(IEL33) = TEMP2(9) + HA(IEL33) 748 HA(IEL33-1) = TEMP2(8) + HA(IEL33-1) 749 HA(IEL33-2) = TEMP2(7) + HA(IEL33-2) 750 HA(IEL33-I*3) = TEMP2(5) + HA(IEL33-I*3) 751 HA(IEL33-I*3-1) = TEMP2(4) + HA(IEL33-I*3-1) 752 HA(IEL33-6*I+1) = TEMP2(1) + HA(IEL33-6*I+1) 753C 754C APPLY SYMMETRY TO DIPOLE TERM AS WELL 755C 756 TEMP(9) = DELDIP(3,K*3 ) 757 TEMP(8) = DELDIP(2,K*3 ) 758 TEMP(7) = DELDIP(1,K*3 ) 759 TEMP(6) = DELDIP(3,K*3-1) 760 TEMP(5) = DELDIP(2,K*3-1) 761 TEMP(4) = DELDIP(1,K*3-1) 762 TEMP(3) = DELDIP(3,K*3-2) 763 TEMP(2) = DELDIP(2,K*3-2) 764 TEMP(1) = DELDIP(1,K*3-2) 765C 766 CALL MAT33(R(1,N), TEMP, TEMP2) 767C 768 DELTMP(3,I*3 ) = TEMP2(9) + DELTMP(3,I*3 ) 769 DELTMP(2,I*3 ) = TEMP2(8) + DELTMP(2,I*3 ) 770 DELTMP(1,I*3 ) = TEMP2(7) + DELTMP(1,I*3 ) 771 DELTMP(3,I*3-1) = TEMP2(6) + DELTMP(3,I*3-1) 772 DELTMP(2,I*3-1) = TEMP2(5) + DELTMP(2,I*3-1) 773 DELTMP(1,I*3-1) = TEMP2(4) + DELTMP(1,I*3-1) 774 DELTMP(3,I*3-2) = TEMP2(3) + DELTMP(3,I*3-2) 775 DELTMP(2,I*3-2) = TEMP2(2) + DELTMP(2,I*3-2) 776 DELTMP(1,I*3-2) = TEMP2(1) + DELTMP(1,I*3-2) 777C 778 30 CONTINUE 779 40 CONTINUE 780C 781 DO 50 I = 1, (NUMAT*3*(NUMAT*3+1))/2 782 50 H(I) = HA(I)/NSYM 783C 784 DO 60 I = 1, 3*NUMAT 785 DELDIP(1,I) = DELTMP(1,I)/NSYM 786 DELDIP(2,I) = DELTMP(2,I)/NSYM 787 60 DELDIP(3,I) = DELTMP(3,I)/NSYM 788C 789 RETURN 790 END 791 SUBROUTINE MAT33(A, B, C) 792C A subroutine that will multiply two 3 by 3 matricies in the following 793C fashion: C = A(transpose) B A 794C 795 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 796 DIMENSION A(9), B(9), C(9), T(9) 797C 798C 799 T(1) = B(1)*A(1) + B(2)*A(4) + B(3)*A(7) 800 T(2) = B(1)*A(2) + B(2)*A(5) + B(3)*A(8) 801 T(3) = B(1)*A(3) + B(2)*A(6) + B(3)*A(9) 802 T(4) = B(4)*A(1) + B(5)*A(4) + B(6)*A(7) 803 T(5) = B(4)*A(2) + B(5)*A(5) + B(6)*A(8) 804 T(6) = B(4)*A(3) + B(5)*A(6) + B(6)*A(9) 805 T(7) = B(7)*A(1) + B(8)*A(4) + B(9)*A(7) 806 T(8) = B(7)*A(2) + B(8)*A(5) + B(9)*A(8) 807 T(9) = B(7)*A(3) + B(8)*A(6) + B(9)*A(9) 808C 809 C(1) = A(1)*T(1) + A(4)*T(4) + A(7)*T(7) 810 C(2) = A(1)*T(2) + A(4)*T(5) + A(7)*T(8) 811 C(3) = A(1)*T(3) + A(4)*T(6) + A(7)*T(9) 812 C(4) = A(2)*T(1) + A(5)*T(4) + A(8)*T(7) 813 C(5) = A(2)*T(2) + A(5)*T(5) + A(8)*T(8) 814 C(6) = A(2)*T(3) + A(5)*T(6) + A(8)*T(9) 815 C(7) = A(3)*T(1) + A(6)*T(4) + A(9)*T(7) 816 C(8) = A(3)*T(2) + A(6)*T(5) + A(9)*T(8) 817 C(9) = A(3)*T(3) + A(6)*T(6) + A(9)*T(9) 818C 819 RETURN 820 END 821 SUBROUTINE SYMP 822 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 823 INCLUDE 'SIZES' 824 COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT 825 COMMON /SYMOPC/ ISYMT(6) 826 CHARACTER*10 ISYMT 827 CHARACTER*5 OPER 828***************************************************************** 829* 830* ON INPUT R = SYMMETRY OPERATIONS (7 MAX) 831* IPO = PERM OPR FOR ABOVE OPERATIONS 832* NSYM = CURRENT NUMBER OF SYMMETRY OPERATIONS 833* NENT = NUMBER OF USER SUPPLIED OPERATIONS 834* 835* ON OUTPUT R = SYMMETRY OPERATIONS (120 MAX) 836* IPO = PERMUTATION OPERATOR FOR SYMMETRY OPERATIONS 837* NSYM = NUMBER OF SYMMETRY OPERATIONS 838* 839***************************************************************** 840C 841C A SUBROUTINE THAT WILL EXPAND THE SYMMETRY OPERATIONS READ IN INTO 842C THE COMPLETE SET. NOTE: VERY FEW OPERATIONS ARE REQUIRED TO 843C GENERATE EVEN VERY LARGE GROUPS OF OPERATIONS. 844C 845C 846 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 847 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 848 2 NCLOSE,NOPEN,NDUMY,FRACT 849 COMMON /COORD / COORD(3,NUMATM) 850 COMMON /KEYWRD/ KEYWRD 851 CHARACTER KEYWRD*241 852C THE NEXT PARAMETERS ARE THE MAX NUMBER OF SYMM FUNCTIONS, THE 853C MAX NUMBER OF SYMM FUNCTIONS TO READ IN, AND THE 854C TOLERENCE USED TO DETERMINE IF TWO FUNCTIONS ARE IDENTICAL 855 PARAMETER (MAXFUN=120) 856 PARAMETER (TOL =1D-2) 857C 858C 859C Variables used: (n represents the number of atomic centers) 860C 861C For the next two items, the last index represents the symmetry 862C operation number. 863C R(9,*): The 9 elements of each record are a packed 3 by 3 864C array of a given symmetry operations. 865C IPO(n,*): A vector that contains the symmetry mapping of atomic ce 866C 867C NSYM IS ALWAYS THE UPPER BOUND OF THE VALID FUNCTIONS. QUIT IF IT 868C REACHES 120. 869C I IS THE SLOW INDEX OF FUNCTIONS TO MULTIPLY 870C J IS THE FAST INDEX OF FUNCTIONS TO MULTIPLY 871C ALWAYS DO R(I)*R(J) AND TAKE I,J FROM 2 TO NSYM 872C 873 I = 2 874 J = 1 875C 876C IF MORE INFORMATION IS WANTED, PRINT HEADDER. 877C 878 IF (INDEX(KEYWRD,' RMAT') .NE. 0)WRITE(6,100) 879 100 FORMAT(/,' ENTERING THE SYMMETRY GENERATING ROUTINE ', 880 +/,' NUMBER SYMM. OPER. * ', 881 + ' NUMBER SYMM. OPER. = ', 882 + ' NUMBER SYMM. OPER.') 883C DETERMINE IF IT IS TIME TO STOP 884C 885 10 J = J+1 886 IF (J .GT. NSYM) THEN 887 J = 2 888 I = I+1 889 IF (I .GT. NSYM) GOTO 50 890 ENDIF 891 IF(NSYM .EQ. MAXFUN) GOTO 50 892C 893C NOW TO START THE MULTIPLICATION 894C 895 R(1,NSYM+1)=R(1,I)*R(1,J)+R(2,I)*R(4,J)+R(3,I)*R(7,J) 896 R(2,NSYM+1)=R(1,I)*R(2,J)+R(2,I)*R(5,J)+R(3,I)*R(8,J) 897 R(3,NSYM+1)=R(1,I)*R(3,J)+R(2,I)*R(6,J)+R(3,I)*R(9,J) 898 R(4,NSYM+1)=R(4,I)*R(1,J)+R(5,I)*R(4,J)+R(6,I)*R(7,J) 899 R(5,NSYM+1)=R(4,I)*R(2,J)+R(5,I)*R(5,J)+R(6,I)*R(8,J) 900 R(6,NSYM+1)=R(4,I)*R(3,J)+R(5,I)*R(6,J)+R(6,I)*R(9,J) 901 R(7,NSYM+1)=R(7,I)*R(1,J)+R(8,I)*R(4,J)+R(9,I)*R(7,J) 902 R(8,NSYM+1)=R(7,I)*R(2,J)+R(8,I)*R(5,J)+R(9,I)*R(8,J) 903 R(9,NSYM+1)=R(7,I)*R(3,J)+R(8,I)*R(6,J)+R(9,I)*R(9,J) 904C 905C IS IT UNIQUE? 906C 907 DO 30 N = 1, NSYM 908 RES = 0.D0 909 DO 20 M = 1, 9 910 20 RES = RES + ABS(R(M,N) - R(M,NSYM+1)) 911 IF (RES .LT. TOL) GOTO 10 912 30 CONTINUE 913C 914C YES, IT IS UNIQUE. NOW, GENERATE THE NEW IPO(,NSYM) 915C 916 NSYM = NSYM + 1 917 DO 40 N = 1, NUMAT 918 40 IPO(N,NSYM) = IPO(IPO(N,J),I) 919C 920C ALL DONE ADDING THE NEW FUNCTION. GO TRY TO FIND A NEW ONE. 921C BUT FIRST, SEE IF WE NEED TO PRINT THIS. 922C 923 IF (INDEX(KEYWRD,' RMAT') .NE. 0) 924 + WRITE(6,110)I,OPER(R(1,I)),J,OPER(R(1,J)),NSYM,OPER(R(1,NSYM)) 925 110 FORMAT(8X,I3,6X,A5,4X,'*',8X,I3,6X,A5,4X,'=',8X,I3,6X,A5) 926 IF (INDEX(KEYWRD,' RMAT') .NE. 0) THEN 927 WRITE(6,120)(R(K,I),K=1,3),(R(K,J),K=1,3),(R(K,NSYM),K=1,3) 928 WRITE(6,130)(R(K,I),K=4,6),(R(K,J),K=4,6),(R(K,NSYM),K=4,6) 929 WRITE(6,140)(R(K,I),K=7,9),(R(K,J),K=7,9),(R(K,NSYM),K=7,9) 930 120 FORMAT(' |',3F7.3,' | |',3F7.3,' | |',3F7.3,' |') 931 130 FORMAT(' |',3F7.3,' | * |',3F7.3,' | = |',3F7.3,' |') 932 140 FORMAT(' |',3F7.3,' | |',3F7.3,' | |',3F7.3,' |',/) 933 ENDIF 934C 935 GOTO 10 936C 937C 938 50 CONTINUE 939C 940C NOW, TO DO FINAL WRAPUP 941C 942 WRITE(6,150)NSYM 943 150 FORMAT(/,' THERE ARE ',I3,' UNIQUE SYMMETRY FUNCTIONS.',/) 944C 945C PRINT THE IPO MATRIX IF ASKED FOR. 946C 947 IF(INDEX(KEYWRD,' IPO') .NE. 0) THEN 948 WRITE(6,160) 949 160 FORMAT(/,20X,'THE PERMUTATION MATRIX') 950 I = 1 951 J = MIN(12,NSYM) 952 60 WRITE(6,170)(K,K=I,J) 953 170 FORMAT(/,/,5X,'OPER. NO. ',12I5) 954 WRITE(6,175)(OPER(R(1,K)),K=I,J) 955 175 FORMAT(5X,'SYMM. OPER. ',12A5) 956 WRITE(6,180) 957 180 FORMAT(5X,'ATOM NO.') 958 DO 70 K = 1, NUMAT 959 70 WRITE(6,190)K,(IPO(K,L),L=I,J) 960 190 FORMAT(I10,5X,12I5) 961 IF (J .LT. NSYM) THEN 962 I = J+1 963 J = MIN(J+12,NSYM) 964 GOTO 60 965 ENDIF 966 ENDIF 967 RETURN 968 END 969 CHARACTER*5 FUNCTION OPER(R) 970C 971 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 972 CHARACTER OPR*5, NUM*10 973 DIMENSION R(9) 974C 975C 976 OPR = ' ' 977 NUM = '0123456789' 978 TRACE = R(1) + R(5) + R(9) 979 DET = R(1)*R(5)*R(9) + R(2)*R(6)*R(7) + R(3)*R(4)*R(8) 980 + - R(1)*R(6)*R(8) - R(2)*R(4)*R(9) - R(3)*R(5)*R(7) 981 TRACE = (TRACE - DET)/2.D0 982 IF (DET .GT. 0.D0) THEN 983 OPR(1:1) = 'C' 984 IF (TRACE .GT. 0.97D0) THEN 985 OPR(1:1) = 'E' 986 GOTO 20 987 ENDIF 988 ELSE 989 OPR(1:1) = 'S' 990 IF (TRACE .GT. 0.97D0) THEN 991 OPR(1:5) = 'Sigma' 992 GOTO 20 993 ENDIF 994 IF (TRACE .LT. -0.97D0) THEN 995 OPR(1:5) = ' Inv ' 996 GOTO 20 997 ENDIF 998 ENDIF 999 IF (TRACE .LT. -0.97D0) THEN 1000 OPR(2:2) = NUM(3:3) 1001 GOTO 20 1002 ENDIF 1003 ANG = ACOS(TRACE) 1004 AFULL = ACOS(-1.0D0)*2.D0 1005 DO 10 I = 3, 18 1006 ANS = I*ANG/AFULL 1007 IF (ABS(ANS - NINT(ANS)) .LE. 2.5D-3) THEN 1008 IF(I .GE.10) THEN 1009 OPR(2:2) = NUM(2:2) 1010 OPR(3:3) = NUM(I-9:I-9) 1011 ELSE 1012 OPR(2:2) = NUM(I+1:I+1) 1013 ENDIF 1014 IF (NINT(ANS) .NE. 1) THEN 1015 OPR(4:5) = '* ' 1016 OPR(5:5) = NUM(NINT(ANS)+1:NINT(ANS)+1) 1017 ENDIF 1018 GOTO 20 1019 ENDIF 1020 10 CONTINUE 1021 OPR(2:5) = 'Unkn' 1022C 1023 20 OPER = OPR 1024 RETURN 1025 END 1026