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