freqcy.f

      SUBROUTINE FREQCY(FMATRX,FREQ,CNORML,REDMAS,TRAVEL,EORC,DELDIP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION FMATRX(*), REDMAS(*), FREQ(*), CNORML(*), TRAVEL(*)
     +,DELDIP(3,*)
      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 /NLLCOM/ FMAT2D(2*MAXPAR**2), VEC(MAXPAR**2)
      COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120)
      COMMON /ATMASS/ ATMASS(NUMATM)
      COMMON /KEYWRD/ KEYWRD
      COMMON /SCRACH/ OLDF(MAXPAR**2)
      CHARACTER KEYWRD*241
      DIMENSION WTMASS(MAXPAR), SHIFT(6), SEC(MAXPAR**2)
      EQUIVALENCE (SEC,OLDF)
      SAVE FACT
      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
      IF(INDEX(KEYWRD,' GROUP').NE.0) THEN
         CALL SYMT(FMATRX, DELDIP)
         IF(INDEX(KEYWRD,' FREQCY').NE.0)THEN
         WRITE(6,'(A)')' SYMMETRIZED HESSIAN MATRIX'
C#         I=-N3
C#         CALL VECPRT(FMATRX,I)
   
C
C   THE FORCE MATRIX IS PRINTED AS AN ATOM-ATOM MATRIX RATHER THAN
C   AS A 3N*3N MATRIX, AS THE 3N MATRIX IS VERY CONFUSING!
C
      IJ=0
      IU=0
      DO 159 I=1,NUMAT
         IL=IU+1
         IU=IL+2
         IM1=I-1
         JU=0
         DO 149 J=1,IM1
            JL=JU+1
            JU=JL+2
            SUM=0.D0
C$DOIT ASIS
            DO 139 II=IL,IU
C$DOIT ASIS
               DO 139 JJ=JL,JU
  139       SUM=SUM+FMATRX((II*(II-1))/2+JJ)**2
            IJ=IJ+1
  149    CNORML(IJ)=SQRT(SUM)
         IJ=IJ+1
  159 CNORML(IJ)=SQRT(
     1FMATRX(((IL+0)*(IL+1))/2)**2+
     2FMATRX(((IL+1)*(IL+2))/2)**2+
     3FMATRX(((IL+2)*(IL+3))/2)**2+2.D0*(
     4FMATRX(((IL+1)*(IL+2))/2-1)**2+
     5FMATRX(((IL+2)*(IL+3))/2-2)**2+
     6FMATRX(((IL+2)*(IL+3))/2-1)**2))
      I=-NUMAT
      CALL VECPRT(CNORML,I)
         ENDIF
      ENDIF
      N3=3*NUMAT
      L=0
      DO 10 I=1,NUMAT
         WEIGHT=1.4142136D0/SQRT(ATMASS(I))
         WTMASS(L+1)=WEIGHT
         WTMASS(L+2)=WEIGHT
         WTMASS(L+3)=WEIGHT
         L=L+3
   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
      I=INDEX(KEYWRD,' K=')
      IF(I.NE.0)THEN
C
C  GO INTO BRILLOUIN ZONE MODE
C
         STEP=READA(KEYWRD,I)
         MONO3=READA(KEYWRD(I:),INDEX(KEYWRD(I:),','))*3
         CALL BRLZON(FMATRX, FMAT2D, N3, SEC, VEC, ALBAND, MONO3, STEP,1
     1)
         RETURN
      ENDIF
       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)-DBLE(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(INDEX(KEYWRD,' GROUP').NE.0) CALL SYMA(FREQ, CNORML)
      IF(EORC) THEN
C
C    CONVERT NORMAL VECTORS TO CARTESIAN COORDINATES
C    (DELETED) AND NORMALIZE SO THAT THE TOTAL MOVEMENT IS 1.0 ANGSTROM.
C
         IJ=0
         DO 110 I=1,N3
            SUM=0.D0
            J=0
            DO 90 JJ=1,NUMAT
               SUM1=0.D0
               CNORML(IJ+1)=CNORML(IJ+1)*WTMASS(J+1)
               SUM1=SUM1+CNORML(IJ+1)**2
C
               CNORML(IJ+2)=CNORML(IJ+2)*WTMASS(J+2)
               SUM1=SUM1+CNORML(IJ+2)**2
C
               CNORML(IJ+3)=CNORML(IJ+3)*WTMASS(J+3)
               SUM1=SUM1+CNORML(IJ+3)**2
C
               J=J+3
               IJ=IJ+3
   90       SUM=SUM+SQRT(SUM1)
            SUM=1.D0/SQRT(SUM)
            IJ=IJ-N3
            DO 100 J=1,N3
               IJ=IJ+1
  100       CNORML(IJ)=CNORML(IJ)*SUM
  110    CONTINUE
C
C          RETURN HESSIAN IN MILLIDYNES/ANGSTROM IN FMATRX
C
         DO 120 I=1,LINEAR
  120    FMATRX(I)=OLDF(I)*1.D-5
      ELSE
C
C  RETURN HESSIAN AS MASS-WEIGHTED FMATRIX
         LINEAR=0
C
         DO 130 I=1,N3
            DO 130 J=1,I
               LINEAR=LINEAR+1
  130    FMATRX(LINEAR)=OLDF(LINEAR)*1.D-5*WTMASS(I)*WTMASS(J)
      ENDIF
      RETURN
      END