SUBROUTINE FREQCY(FMATRX,FREQ,CNORML,REDMAS,TRAVEL,EORC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE "SIZES"
      DIMENSION FMATRX(*), REDMAS(*), FREQ(*), CNORML(*), TRAVEL(*)
      LOGICAL EORC
*********************************************************************
*
*  FRCE CALCULATES THE FORCE CONSTANTS AND VIBRATIONAL FREQUENCIES
*       FOR A MOLECULE.  IT USES THE ISOTOPIC MASSES TO WEIGHT THE
*       FORCE MATRIX
*
* ON INPUT   FMATRX   =  FORCE MATRIX, OF SIZE NUMAT*3*(NUMAT*3+1)/2.
*
*********************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /ATMASS/ ATMASS(NUMATM)
      COMMON /SCRACH/ OLDF(MAXHES), DUMMY(MAXPAR**2-MAXHES)
      DIMENSION WTMASS(MAXPAR), SHIFT(6)
      DATA FACT/6.023D23/
C
C    CONVERSION FACTOR FOR SPEED OF LIGHT AND 2 PI.
C
      C2PI=1.D0/(2.998D10*3.141592653598D0*2.D0)
C NOW TO CALCULATE THE VIBRATIONAL FREQUENCIES
C
C   FIND CONVERSION CONSTANTS FOR MASS WEIGHTED SYSTEM
      N3=3*NUMAT
      L=0
      DO 10 I=1,NUMAT
         WEIGHT=1.4142136D0/SQRT(ATMASS(I))
         DO 10 J=1,3
            L=L+1
   10 WTMASS(L)=WEIGHT
C    CONVERT TO MASS WEIGHTED FMATRX
      LINEAR=0
      DO 20 I=1,N3
         DO 20 J=1,I
            LINEAR=LINEAR+1
            OLDF(LINEAR)=  FMATRX(LINEAR)*1.D5
   20 FMATRX(LINEAR)=FMATRX(LINEAR)*WTMASS(I)*WTMASS(J)
C
C    1.D5 IS TO CONVERT FROM MILLIDYNES/ANGSTROM TO DYNES/CM.
C
C    DIAGONALIZE
      CALL FRAME(FMATRX,NUMAT,1, SHIFT)
      CALL RSP(FMATRX,N3,N3,FREQ,CNORML)
      DO 30 I=1,N3
         J=(FREQ(I)+50.D0)*0.01D0
   30 FREQ(I)=FREQ(I)-J*100
      DO 40 I=1,N3
   40 FREQ(I)=FREQ(I)*1.D5
C
C    CALCULATE REDUCED MASSES, STORE IN REDMAS
C
      DO 80 I=1,N3
         II=(I-1)*N3
         SUM=0.D0
         DO 70 J=1,N3
            JII=J+II
            JJ=(J*(J-1))/2
            DO 50 K=1,J
   50       SUM=SUM+CNORML(JII)*OLDF(JJ+K)*CNORML(K+II)
            DO 60 K=J+1,N3
   60       SUM=SUM+CNORML(JII)*OLDF((K*(K-1))/2+J)*CNORML(K+II)
   70    CONTINUE
         SUM1=SUM*2.D0
         IF(ABS(FREQ(I)).GT.ABS(SUM)*1.D-20) THEN
            SUM=1.D0*SUM/FREQ(I)
         ELSE
            SUM=0.D0
         ENDIF
         FREQ(I)=SIGN(SQRT(FACT*ABS(FREQ(I)))*C2PI,FREQ(I))
         IF(ABS(FREQ(I)).LT.ABS(SUM1)*1.D+20) THEN
            SUM1=SQRT(ABS(FREQ(I)/(SUM1*1.D-5)))
         ELSE
            SUM1=0.D0
         ENDIF
         IF(SUM.LT.0.D0.OR.SUM.GT.100)SUM=0.D0
C
C 0.0063024=SQRT(2*A*B*C/N) WHERE
C         A=1.196D8 = CONVERSION OF CM**(-1) TO (ERGS = DYNE.ANGSTROMS)
C         B=1000.0  = MILLIDYNES TO DYNES
C         C=1.D8    = CENTIMETERS TO ANGSTROMS
C         N=6.02205D23 = AVOGADRO'S NUMBER
         TRAVEL(I)=SUM1*0.0063024D0
         IF(TRAVEL(I).GT.1.D0)TRAVEL(I)=0.D0
C#      WRITE(6,*)TRAVEL(I)
   80 REDMAS(I)=SUM
      IF(EORC) THEN
C
C    SWITCH EIGENVALUES TO FREQUENCIES
C
C    CONVERT NORMAL VECTORS TO CARTESIAN COORDINATES
C    AND NORMALIZE SO THAT THE TOTAL MOVEMENT IS 1.0 ANGSTROM.
C
         IJ=0
         DO 120 I=1,N3
            SUM=0.D0
            J=0
            DO 100 JJ=1,NUMAT
               SUM1=0.D0
               DO 90 J1=1,3
                  J=J+1
                  IJ=IJ+1
                  CNORML(IJ)=CNORML(IJ)*WTMASS(J)
   90          SUM1=SUM1+CNORML(IJ)**2
  100       SUM=SUM+SQRT(SUM1)
            SUM=1.D0/SUM
            IJ=IJ-N3
            DO 110 J=1,N3
               IJ=IJ+1
  110       CNORML(IJ)=CNORML(IJ)*SUM
  120    CONTINUE
C
C          RETURN HESSIAN IN MILLIDYNES/ANGSTROM IN FMATRX
C
         DO 130 I=1,LINEAR
  130    FMATRX(I)=OLDF(I)*1.D-5
      ELSE
C
C  RETURN HESSIAN AS MASS-WEIGHTED FMATRIX
         LINEAR=0
C
         DO 140 I=1,N3
            DO 140 J=1,I
               LINEAR=LINEAR+1
  140    FMATRX(LINEAR)=OLDF(LINEAR)*1.D-5*WTMASS(I)*WTMASS(J)
      ENDIF
      RETURN
      END