C PROGRAM GEOMET(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) C PURPOSE: (1) TO READ IN BOND LENGTHS AND ANGLES AND OUTPUT CARTESIAN C COORDINATES. C (2) TO CALCULATE MOMENTS OF INERTIA C (3) TO CALCULATE HINDRANCES. C C CALCULATION OF CARTESIAN COORDINATES USES METHOD OF H.B. THOMPSON, C J. CHEM. PHYS. 47, 3407 (1967). C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL COORRD,RDRAD DOUBLE PRECISION MASS(50) DIMENSION ITITLE(20),NA(50),NB(50),A(50),B(50),R(50),LISTA(50), 1 AJ(50,16),X(50),Y(50),Z(50),DIST(50),SMASS(50),AORD(3) 2,AEIG(9),REIG(9),AM(7),XS(50),YS(50),ZS(50),COORD(50,3), 3 ATVR(50),ATR(50),LATR(50) COMMON /TRANS/ BB(16),F DATA EMX/0./,EMY/0./,EMZ/0./,EM/0./,AEIG/9*0./,ATR/50*.0/ TWOPI=8.*ATAN(1.) F=TWOPI/360. BB(4)=0. BB(8)=0. BB(9)=0. BB(12)=0. BB(16)=1. READ(5,400) ITITLE WRITE(6,401) ITITLE 400 FORMAT(20A4) 401 FORMAT(/,20X,20A4,/,10X,'LATEST UPDATE JULY 25 1991') READ(5,*) NN COORRD = NN.LE.0 NN=IABS(NN) WRITE(6,1) NN 1 FORMAT(' NUMBER OF ATOMS:',I5) NATOMS=NN IF(COORRD) GO TO 872 WRITE(6,404) 404 FORMAT(/,' ATOM NO. ATTACHED TO BOND DIHEDRAL', 1 T43,'BOND LENGTH',T58,'MASS',/ 1 ' ANGLE ANGLE', 2 /,T23,'(DEGREES) (DEGREES)',T44,'(ANGSTROM)',T58,'(AMU)') NN=NN+1 NA(1)=1 NB(1)=0 A(1)=0. B(1)=0. R(1)=0. DO 15 I=2,NN READ(5,*) NA(I),NB(I),A(I),B(I),R(I),MASS(I) IF(I.EQ.2.OR.I.EQ.3) A(I)=0 IF(I.EQ.2.OR.I.EQ.3) B(I)=0 IF(I.EQ.2) R(I)=0. NAI=NA(I) NA(I)=NAI+1 NBI=NB(I) NB(I)=NBI+1 WRITE(6,405) NAI,NBI,A(I),B(I),R(I),MASS(I) 405 FORMAT(I5,I10,5X,2F10.2,T41,F10.3,T55,F8.3) 15 CONTINUE DO 20 I=1,16 AJ(1,I)=0. 20 CONTINUE DO 36 I=2,NN AI=A(I) BI=B(I) RI=R(I) CALL LOADB(AI,BI,RI) IF(NB(I)-1) 50,40,50 40 DO 41 II=1,16 AJ(I,II)=BB(II) 41 CONTINUE GO TO 36 50 N=NB(I) DO 151 II=1,4 DO 51 NJ=1,4 SUM=0. DO 52 JJ=1,4 SUM=SUM+AJ(N,II+(JJ-1)*4)*BB((NJ-1)*4+JJ) 52 CONTINUE AJ(I,II+4*(NJ-1))=SUM 51 CONTINUE 151 CONTINUE 36 CONTINUE GO TO 873 872 CONTINUE NN=NN+1 DO 870 II=2,NN READ(5,*) AJ(II,13),AJ(II,14),AJ(II,15),MASS(II) NA(II)=II 870 CONTINUE 873 CONTINUE WRITE(6,490) 490 FORMAT(//,10X,'CARTESIAN COORDINATES:',//, 1 ' ATOM MASS X Y Z',/) DO 70 II=2,NN X(II)=AJ(II,13) Y(II)=AJ(II,14) Z(II)=AJ(II,15) AMM=MASS(II) EM=EM+AMM EMX=EMX+AMM*X(II) EMY=EMY+AMM*Y(II) EMZ=EMZ+AMM*Z(II) WRITE(6,500) II-1,MASS(II),X(II),Y(II),Z(II) 500 FORMAT(I5,F10.3,3F10.4) 70 CONTINUE EMX=EMX/EM EMY=EMY/EM EMZ=EMZ/EM WRITE(6,779) (I-1,I=2,NN) 779 FORMAT(//,10X,'DISTANCES:',//,2X,4(15I5,/)) DO 773 II=3,NN XX=X(II) YY=Y(II) ZZ=Z(II) DO 778 J=2,II-1 XLOC=X(J) YLOC=Y(J) ZLOC=Z(J) DIST(J)=DSQRT((XX-XLOC)**2+(YY-YLOC)**2+(ZZ-ZLOC)**2) IF(DIST(J).GT.0.95) GO TO 778 WRITE(6,76) II-1,J-1 76 FORMAT(' WARNING, ATOMS ',I3,' AND ',I3,' ARE VERY CLOSE!') 778 CONTINUE WRITE(6,772) II-1,(DIST(J),J=2,II-1) 772 FORMAT(I3,4(15F5.2,/)) 773 CONTINUE DO 3 II=2,NN XI=X(II)-EMX YI=Y(II)-EMY ZI=Z(II)-EMZ AMM=MASS(II) RI2=XI*XI + YI*YI + ZI*ZI AEIG(1)=AEIG(1) + AMM*(RI2-XI*XI) AEIG(2)=AEIG(2) - AMM*XI*YI AEIG(3)=AEIG(3) + AMM*(RI2-YI*YI) AEIG(4)=AEIG(4) - AMM*XI*ZI AEIG(5)=AEIG(5) - AMM*YI*ZI AEIG(6)=AEIG(6) + AMM*(RI2-ZI*ZI) 3 CONTINUE ITHREE=3 NOUGHT=0 CALL EIGEN(AEIG,REIG,ITHREE,NOUGHT) C C PUT MOMENTS OF INERTIA IN ORDER SUCH THAT IC C IS THE 'MOST UNEQUAL' (SIC) ONE C AORD(1)=AEIG(1) AORD(2)=AEIG(3) AORD(3)=AEIG(6) DO 142 I=1,3 DO 43 J=I,3 I3=3*I I2=I3-1 I1=I3-2 J3=3*J J2=J3-1 J1=J3-2 IF(AORD(I).LE.AORD(J)) GO TO 43 TEMPA=AORD(I) TEMPR1=REIG(I1) TEMPR2=REIG(I2) TEMPR3=REIG(I3) AORD(I)=AORD(J) REIG(I1)=REIG(J1) REIG(I2)=REIG(J2) REIG(I3)=REIG(J3) AORD(J)=TEMPA REIG(J1)=TEMPR1 REIG(J2)=TEMPR2 REIG(J3)=TEMPR3 43 CONTINUE 142 CONTINUE R21=AORD(2)/AORD(1) R32=AORD(3)/AORD(2) IF(R32.GE.R21) GO TO 149 TEMPA=AORD(3) TEMPR1=REIG(7) TEMPR2=REIG(8) TEMPR3=REIG(9) AORD(3)=AORD(1) REIG(7)=REIG(1) REIG(8)=REIG(2) REIG(9)=REIG(3) AORD(1)=TEMPA REIG(1)=TEMPR1 REIG(2)=TEMPR2 REIG(3)=TEMPR3 149 CONTINUE C C CALCULATE CARTESIAN COORDINATES IN PRINCIPAL AXES SYSTEM C WRITE(6,345) 345 FORMAT(//,2X,'CARTESIAN COORDINATES IN PRINCIPAL AXES SYSTEM: ' *,//,' ATOM MASS X Y Z',/) DO 27 II=2,NN XI=X(II)-EMX YI=Y(II)-EMY ZI=Z(II)-EMZ XPRINC=XI*REIG(1)+YI*REIG(2)+ZI*REIG(3) YPRINC=XI*REIG(4)+YI*REIG(5)+ZI*REIG(6) ZPRINC=XI*REIG(7)+YI*REIG(8)+ZI*REIG(9) WRITE(6,500) II-1,MASS(II),XPRINC,YPRINC,ZPRINC 27 CONTINUE WRITE(6,28) AEIG(1),AEIG(3),AEIG(6) 28 FORMAT(//,' PRINCIPAL MOMENTS OF INERTIA (AMU ANGSTROM**2)',/, 1 3F12.3) AMM=(AEIG(1)*AEIG(3)*AEIG(6))**0.333333333 WRITE(6,29) AMM 29 FORMAT(/,' IM = (IA*IB*IC)**(1/3)=',F12.3) READ(5,*) ITEST IF(ITEST.EQ.0) STOP WRITE(6,23) 23 FORMAT(//,20X,'CALCULATION OF HINDRANCES') READ(5,*) NATOMA,INDEXA RDRAD = NATOMA.LT.0 NATOMA = IABS(NATOMA) READ(5,*) (LISTA(I),I=1,NATOMA) READ(5,*) INDEXB WRITE(6,56) NATOMA,INDEXA,INDEXB 56 FORMAT(/,' NUMBER OF ATOMS IN MOIETY A =',T37,I5,/, 1 ' INDEX OF PIVOT ATOM IN A =',T37,I5,/, 2 ' INDEX OF PIVOT ATOM IN B =',T37,I5,/) IF(.NOT.RDRAD) GO TO 8341 READ(5,*) IVAN WRITE (6,53) 53 FORMAT(/,' INDEXES AND RADII OF ATOMS SELECTED FOR RADII INPUT ') DO 127 I= 1,IVAN READ (5,*) LATR(I),ATVR(I) ATR(LATR(I))= ATVR(I) WRITE(6,54) LATR(I),ATVR(I) 54 FORMAT(10X,I4,F12.3) WRITE(6,*) 127 CONTINUE 8341 CONTINUE WRITE(6,97) (LISTA(I),MASS(LISTA(I)+1),I=1,NATOMA) 97 FORMAT(' INDEXES OF ATOMS IN MOIETY A:',T40,'INDEX',T50,'MASS', 1 /,50(T40,I3,T47,F6.1,/)) C STORE COORDINATES IN XS, YS AND ZS SO THAT ATOM 1 IS PIVOT FOR FIRST C MOIETY, ATOM (NATOMA+1) IS PIVOT FOR SECOND MOIETY, AND OTHER ATOMS OF C FIRST MOIETY ARE ARRANGED 2 ... NATOMA, OTHER ATOMS OF SECOND ARE ARRANGED C NATOMA+1 ... NATOMS C XOLD=X(INDEXA+1) YOLD=Y(INDEXA+1) ZOLD=Z(INDEXA+1) XS(1)=X(INDEXA+1) YS(1)=Y(INDEXA+1) ZS(1)=Z(INDEXA+1) SMASS(1)=MASS(INDEXA+1) XS(NATOMA+1)=X(INDEXB+1) YS(NATOMA+1)=Y(INDEXB+1) ZS(NATOMA+1)=Z(INDEXB+1) SMASS(NATOMA+1)=MASS(INDEXB+1) IC=1 DO 326 I=1,NATOMA IF(LISTA(I).EQ.INDEXA) GO TO 326 IC=IC+1 IA=LISTA(I) XS(IC)=X(IA+1) YS(IC)=Y(IA+1) ZS(IC)=Z(IA+1) SMASS(IC)=MASS(IA+1) 326 CONTINUE IC=NATOMA+1 DO 327 II=2,NN I=II-1 IF(I.EQ.INDEXA.OR.I.EQ.INDEXB) GO TO 327 DO 328 J=1,NATOMA IF(I.EQ.LISTA(J)) GO TO 327 328 CONTINUE IC=IC+1 XS(IC)=X(II) YS(IC)=Y(II) ZS(IC)=Z(II) SMASS(IC)=MASS(II) 327 CONTINUE C TRANSFORM SO THAT PIVOT OF A IS ORIGIN DO 30 I=1,NATOMS X(I)=XS(I)-XOLD Y(I)=YS(I)-YOLD Z(I)=ZS(I)-ZOLD 30 CONTINUE C NOW ROTATE ABOUT XZ, THEN ABOUT YZ, SO THAT PIVOTS ARE ON AXIS. XA=X(NATOMA+1) YA=Y(NATOMA+1) ZA=Z(NATOMA+1) THETA=ATAN2(-XA,ZA) C=DCOS(THETA) S=DSIN(THETA) DO 32 I=2,NATOMS XOLD=X(I) ZOLD=Z(I) X(I)=C*XOLD + S*ZOLD Z(I)=-S*XOLD + C*ZOLD 32 CONTINUE ZA=Z(NATOMA+1) THETA=ATAN2(-YA,ZA) C=DCOS(THETA) S=DSIN(THETA) DO 33 I=2,NATOMS YOLD=Y(I) ZOLD=Z(I) Y(I)=C*YOLD + S*ZOLD Z(I)=-S*YOLD + C*ZOLD 33 CONTINUE DO 67 I=1,NATOMS COORD(I,1)=X(I) COORD(I,2)=Y(I) COORD(I,3)=Z(I) 67 CONTINUE CALL HINDRT(TWOPI,NATOMS,NATOMA,SMASS,COORD,ATR) STOP END SUBROUTINE LOADB(AA,BI,RR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /TRANS/ BB(16),F AB=F*AA BC=F*BI CA=DCOS(AB) SA=DSIN(AB) SB=DSIN(BC) CB=DCOS(BC) BB(1)=-CA BB(2)=SA*CB BB(3)=SA*SB BB(5)=-SA BB(6)=-CA*CB BB(7)=-CA*SB BB(10)=-SB BB(11)=CB BB(13)=-RR*CA BB(14)=RR*SA*CB BB(15)=RR*SA*SB RETURN END SUBROUTINE HINDRT(TWOPI,NATOMS,NATOMA,MASS,COORD,ATR) C FINDS HINDRANCE ANGLE FOR TWO MOIETIES, TOGETHER WITH MOMENTS OF INERTIA. C AUTHORS: P G GREENHILL, S C SMITH, R G GILBERT, A R WHYTE, C M J JORDAN, I G PITT IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION MASS(50),COORD(50,3),NCOORD(50,3),RAD(50), 1 INCTHETA,INCPHI,AVHIND,THETA,PHI,HIND(400),TEMPX(50),TEMPY(50), 2 TEMPZ(50),R(50),DIST,MIN,AEIG(9),COORD2(50,3),MASS2(50),COFM(3), 3 AORD(3),THIND,ATR(50),RAD2(50) INTEGER NATOMS,NATOMA,ATOM,AXIS,B,A,M,CX,Q,PATOM LOGICAL FLAG,NOTORS(2),NO2D(2) DATA NANG/362/ C- C- NATOMS=TOTAL NUMBER OF ATOMS ; NATOMA=NUMBER OF ATOMS IN MOIETY C FIRST ATOM IS PIVOT. C- FEED IN COORDS FOR ATOMS ON AXIS AT 1 AND NATOMA+1 TORS=0. DO 224 I224=1,2 NOTORS(I224)=.FALSE. NO2D(I224)=.FALSE. 224 CONTINUE II5=0 WRITE(6,9985) 9985 FORMAT(/,' ATOM MASS RADIUS') DO 9 ATOM=1,NATOMS RAD(ATOM) = 0. M=INT(MASS(ATOM)+0.5) C ASSIGN VAN DER WAALS RADIUS FROM ATOMIC WEIGHT IF (M .EQ. 12) RAD(ATOM)=1.72 IF (M.LE. 3) RAD(ATOM)=1.27 IF (M .EQ. 16) RAD(ATOM)=1.69 IF (M .EQ. 14) RAD(ATOM)=1.69 IF (M .EQ. 19) RAD(ATOM)=1.67 IF (M.GE.35.AND.M.LE.38) RAD(ATOM)=1.94 IF(M.EQ.32) RAD(ATOM)=1.99 IF(M.EQ.80) RAD(ATOM)=2.09 IF(M.EQ.127) RAD(ATOM)=2.28 C FOLLOWING ARE (IN ORDER) P, AS, SI AND GE, OBTAINED BY C BENSON'S METHOD OF ADDING 0.95 A TO COVALENT RADIUS. IF(M.EQ.31) RAD(ATOM)=2.05 IF(M.EQ.75) RAD(ATOM)=2.17 IF(M.EQ.28) RAD(ATOM)=2.12 IF(M.EQ.73) RAD(ATOM)=2.16 C NEXT ARE SB, SE IF(M.EQ.122) RAD(ATOM)=2.36 IF(M.EQ.79) RAD(ATOM)=2.12 IF(RAD(ATOM).NE.0. .OR. ATR(ATOM).NE.0.) GO TO 9988 WRITE (6,9987) 9987 FORMAT (/,' NO VDW RADIUS IN LIST FOR THIS ATOMIC MASS',/, ' USE $ INPUT OF RADIUS OPTION FOR THIS ATOM') STOP 9988 IF(ATR(ATOM).GT.0.) RAD(ATOM) = ATR(ATOM) WRITE(6,9986) ATOM,MASS(ATOM),RAD(ATOM) 9986 FORMAT(I3,2F10.3) 9 CONTINUE DO 555 PATOM=1,2 IF(PATOM.EQ.1)GOTO 521 DO 503 Q=NATOMA+1,NATOMS CX=Q-NATOMA DO 504 AXIS=1,3 COORD2(CX,AXIS)=COORD(Q,AXIS) 504 CONTINUE MASS2(CX)=MASS(Q) RAD2(CX)=RAD(Q) 503 CONTINUE DO 505 Q=1,NATOMA CX=(NATOMS-NATOMA)+Q DO 506 AXIS=1,3 COORD2(CX,AXIS)=COORD(Q,AXIS) 506 CONTINUE MASS2(CX)=MASS(Q) RAD2(CX)=RAD(Q) 505 CONTINUE DO 507 Q=1,NATOMS COORD(Q,1)=COORD2(Q,1) COORD(Q,2)=-COORD2(Q,2) COORD(Q,3)=-COORD2(Q,3) MASS(Q)=MASS2(Q) RAD(Q)=RAD2(Q) 507 CONTINUE NATOMA=NATOMS-NATOMA WRITE(6,508) 508 FORMAT(//,T10,'***** ROTATION ABOUT OTHER PIVOT ATOM *****',//) 521 CONTINUE WRITE(6,89) C- C- NCOORDS IS ARRAY WITH ATOMIC COORDS RELATIVE TO PIVOT ATOM C- DO 10 ATOM=1,NATOMS DO 11 AXIS=1,3 NCOORD(ATOM,AXIS)=COORD(ATOM,AXIS)-COORD(1,AXIS) 11 CONTINUE 10 CONTINUE C- C- INCTHETA=1 DEGREE, INCPHI=10 DEGREES. C- INCTHETA=TWOPI/360. INCPHI=10.*INCTHETA C- C- THE FOLLOWING FINDS HINDRANCE IN THE THETA DIRECTION AS A C- FUNCTION OF PHI - WHERE THETA AND PHI ARE THE EULER ANGLES C- DESCRIBING FRAGMENT B. C- TOTAL HINDRANCE IS THEN AN AVERAGE OVER PHI. C- NOTE THAT THE FOLLOWING ALGORITHM DOES NOT AVERAGE OVER ALL C- POSSIBLE ORIENTATIONS OF THE TWO FRAGMENTS, BUT ONLY OVER C- THEIR RELATIVE ORIENTATION DESCRIBED BY THE EULER ANGLE PHI. C- THE NEGLECTED CONFIGURATIONS ARE THOSE ROTATIONS THAT OCCUR C- WITH THE ORIENTATION OF THE TWO FRAGMENTS HELD CONSTANT. C- AVHIND=0. PHI=0. DO 777 IPHI=1,351,10 HIND(IPHI)=0. THETA=0. DO 700 ITHETA=1,181 C- C- TEMPX ETC ARE THE TEMPORARY COORDINATE VALUES FORMED BY THE C- ROTATIONS IN PHI AND THETA. C- DO 13 B=NATOMA+1,NATOMS TEMPX(B)=NCOORD(B,1)*DCOS(PHI)+NCOORD(B,2)*DSIN(PHI) TEMPY(B)=-NCOORD(B,1)*DCOS(THETA)*DSIN(PHI)+NCOORD(B,2)* 1 DCOS(THETA)*DCOS(PHI)+NCOORD(B,3)*DSIN(THETA) TEMPZ(B)=NCOORD(B,1)*DSIN(THETA)*DSIN(PHI)-NCOORD(B,2)* 1 DSIN(THETA)*DCOS(PHI)+NCOORD(B,3)*DCOS(THETA) 13 CONTINUE C- C- CHECK TO SEE IF ATOMS ARE CLOSER THAN V DER W RADII C- DO 14 A=1,NATOMA DO 15 B=NATOMA+1,NATOMS IF ((A.EQ.1) .AND. (B.EQ.NATOMA+1)) GOTO 15 DIST=DSQRT((TEMPX(B)-NCOORD(A,1))**2 + 1 (TEMPY(B)-NCOORD(A,2))**2 + (TEMPZ(B)-NCOORD(A,3))**2) MIN=RAD(A)+RAD(B) IF (DIST.LE.MIN) GOTO 111 15 CONTINUE 14 CONTINUE THETA=THETA+INCTHETA 700 CONTINUE 111 HIND(IPHI)=THETA IF((HIND(IPHI)).GT.(0.5*TWOPI))HIND(IPHI)=0.5*TWOPI HIND(IPHI)=HIND(IPHI)/(TWOPI/INCPHI) AVHIND=AVHIND+HIND(IPHI) PHI=PHI+INCPHI 777 CONTINUE 89 FORMAT(/,20X,20('*'),/) C- IDEG=INT((AVHIND*360./TWOPI)+0.5) IF(IDEG.LE.179) GO TO 62 IDEG=180 AVHIND=0.5*TWOPI WRITE(6,61) 61 FORMAT(' MOIETIES ARE UNHINDERED') GO TO 63 62 CONTINUE THIND=2.0*AVHIND ITDEG=2*IDEG WRITE(6,90) THIND,ITDEG 90 FORMAT(/,' AS AN AVERAGE OVER THE EULER COORDINATE PHI',/, 1 ' TOTAL HINDRANCE ANGLE =',F10.3,' RADIANS =',I4,' DEGREES') WRITE(6,91)AVHIND,IDEG 91 FORMAT(/,' THE AVERAGE HINDRANCE (OVER PHI) IN THE EULER ',/, 1 ' COORDINATE THETA IS ',F10.3,' RADIANS =',I4,' DEGREES') 63 CONTINUE WRITE(6,89) C- C- GENERATE MOMENTS OF INERTIA ABOUT CENTER OF MASS OF MOIETY C- C- FIND CENTER OF MASS OF MOIETY C- DO 104 I=1,3 COFM(I)=0. 104 CONTINUE DO 339 Q=1,9 AEIG(Q)=0.0 339 CONTINUE C FIND MOMENTS OF INERTIA ON AXIS OF PIVOT DO 3 I=1,NATOMA XI=NCOORD(I,1) YI=NCOORD(I,2) ZI=NCOORD(I,3) AMM=MASS(I) RI2=XI*XI + YI*YI + ZI*ZI AEIG(1)=AEIG(1) + AMM*(RI2-XI*XI) AEIG(2)=AEIG(2) - AMM*XI*YI AEIG(3)=AEIG(3) + AMM*(RI2-YI*YI) AEIG(4)=AEIG(4) - AMM*XI*ZI AEIG(5)=AEIG(5) - AMM*YI*ZI AEIG(6)=AEIG(6) + AMM*(RI2-ZI*ZI) 3 CONTINUE ITHREE=3 NOUGHT=0 CALL EIGEN(AEIG,R,ITHREE,NOUGHT) WRITE(6,28) AEIG(1),AEIG(3),AEIG(6) 28 FORMAT(//,' PRINCIPAL MOMENTS OF INERTIA AROUND ', 1 'PIVOT', 2 /,3F12.3,' (AMU ANGSTROM**2)') C PUT MOMENTS OF INERTIA IN ORDER AORD(1)=AEIG(1) AORD(2)=AEIG(3) AORD(3)=AEIG(6) DO 142 I=1,3 DO 43 J=I,3 IF(AORD(I).LE.AORD(J)) GO TO 43 TEMP=AORD(I) AORD(I)=AORD(J) AORD(J)=TEMP 43 CONTINUE 142 CONTINUE C PUT IN ORDER SO THAT LAST MOMENT OF INERTIA IS THE UNEQUAL ONE (IC) DO 921 I21=1,3 IF(AORD(I21).LT..05)AORD(I21)=.05 921 CONTINUE R21=AORD(2)/AORD(1) R32=AORD(3)/AORD(2) IF(R32.GE.R21) GO TO 149 TEMP=AORD(3) AORD(3)=AORD(1) AORD(1)=TEMP 149 CONTINUE IF(AORD(3).EQ..05)NOTORS(PATOM)=.TRUE. IF(NOTORS(PATOM))II5=PATOM IF((AORD(1).EQ..05).OR.(AORD(2).EQ..05))NO2D(PATOM)=.TRUE. TORS=TORS+(1./AORD(3)) BVAL=DSQRT(AORD(1)*AORD(2)) IF (AVHIND.EQ.0.0) GO TO 230 IF(AVHIND.LE.(0.5*TWOPI))BHIND=BVAL*0.5*(1.-DCOS(AVHIND)) BVAL=16.85/BVAL IF((AVHIND.LE.(0.5*TWOPI)))BHIND=16.85/BHIND IF((.NOT.NO2D(PATOM)).AND.(AVHIND.LT.(0.5*TWOPI))) 1 WRITE(6,276) BHIND 276 FORMAT(/,15X,43('-'),/,10X,'ROTATIONAL CONSTANT B FOR 2-DIMENS', 1 'IONAL HINDERED ROTOR =',/,10X,1P,E10.3,' CM-1', 2 ' (COS THETA HINDRANCE TERM INCLUDED)',/,15X,43('-'),/) GO TO 277 230 BVAL=16.85/BVAL WRITE(6,231) 231 FORMAT(/,' HINDERED ROTOR TREATMENT IS INAPPLICABLE',/,' TO', 1 ' ROTATIONS WITH HINDRANCE ANGLE ZERO',/) 277 IF(.NOT.NO2D(PATOM))WRITE(6,278) BVAL 278 FORMAT(/,15X,43('-'),/,14X,'ROTATIONAL CONSTANT B FOR 2-DIMENS', 1 'IONAL FREE ROTOR =',/,10X,1P,E10.3,' CM-1',/,15X,43('-'),/) IF(NO2D(PATOM))WRITE(6,222)PATOM 222 FORMAT(/,1X,'NO TWO-DIMENSIONAL ROTOR FOR MOIETY ',I2, 1 ' : MOMENT OF INERTIA IS ZERO.',/) 555 CONTINUE BVAL=16.85*TORS IF((.NOT.NOTORS(1)).AND.(.NOT.NOTORS(2)))WRITE(6,432)BVAL 432 FORMAT(//,10X,' ROTATIONAL CONSTANT B FOR OVERALL TORSION =', 1 1P,E10.3,' CM-1') IF(NOTORS(1).OR.NOTORS(2))WRITE(6,223)II5 223 FORMAT(/,1X,'NO OVERALL TORSIONAL DEGREE OF FREEDOM : ', 1 'MOMENT OF INERTIA OF MOIETY ',I2,' ABOUT THE UNIQUE', 2 ' AXIS IS ZERO.',/) RETURN END SUBROUTINE EIGEN(A,R,N,MV) C DESCRIPTION OF PARAMETERS C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF C MATRIX A IN DESCENDING ORDER. IMPLICIT DOUBLE PRECISION (A-H,O-Z) C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, C IN SAME SEQUENCE AS EIGENVALUES) C N - ORDER OF MATRICES A AND R C MV- INPUT CODE C 0 COMPUTE EIGENVALUES AND EIGENVECTORS C 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE C DIMENSIONED BUT MUST STILL APPEAR IN CALLING C SEQUENCE) C DIMENSION A(1),R(1) C C ............................................................... C C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION C STATEMENT WHICH FOLLOWS. C C DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, C 1 COSX2,SINCS,RANGE C C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. DSQRT IN STATEMENTS C 40, 68, 75, AND 78 MUST BE CHANGED TO DDSQRT. DABS IN STATEMENT C 62 MUST BE CHANGED TO DDABS. THE CONSTANT IN STATEMENT 5 SHOULD C BE CHANGED TO 1.0D-12. C C ............................................................... C C GENERATE IDENTITY MATRIX C 5 RANGE=1.0E-6 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0 IF(I-J) 20,15,20 15 R(IJ)=1.0 20 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) C 25 ANORM=0.0 DO 35 I=1,N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM) 165,165,40 40 ANORM=1.414*DSQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR C IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 C C COMPUTE SIN AND COS C 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF( DABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=0.5*(A(LL)-A(MM)) 68 Y=-A(LM)/ DSQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/ DSQRT(2.0*(1.0+( DSQRT(1.0-Y*Y)))) SINX2=SINX*SINX 78 COSX= DSQRT(1.0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLETION C C TEST FOR M = LAST COLUMN C 130 IF(M-N) 135,140,135 135 M=M+1 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C 160 IF(THR-ANRMX) 165,165,45 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE RETURN END