      SUBROUTINE FORCE
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      PARAMETER (MAXPA2=MAXPAR**2)
      common /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY,DUMY(MAXPAR)
      common /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),
     1                     LOCDEP(MAXPAR)
      common /GEOKST/ NATOMS,LABELS(NUMATM),
     1NA(NUMATM),NB(NUMATM),NC(NUMATM)
      common /FMATRX/ FMATRX(MAXHES)
      common /KEYWRD/ KEYWRD
      common /GRADNT/ GRAD(MAXPAR),GNORM
      PARAMETER (IPADD=2*MORB2+2*MAXORB-MAXPAR-MAXPA2)
      common /VECTOR/ CNORML(MAXPA2),FREQ(MAXPAR),DUMMY(IPADD)
      common /ELEMTS/ ELEMNT(107)
      common /LAST  / LAST
      common /MESAGE/ IFLEPO,ISCF
      common /SIMBOL/ SIMBOL(MAXPAR)
      common /GEOM  / GEO(3,NUMATM)
      common /COORD / COORD(3,NUMATM)
***********************************************************************
*
*   FORCE CALCULATES THE FORCE CONSTANTS FOR THE MOLECULE, AND THE
*         VIBRATIONAL FREQUENCIES.  ISOTOPIC SUBSTITUTION IS ALLOWED.
*
***********************************************************************
      common /EULER / TVEC(3,3), ID
      common /SCFTYP/ EMIN, LIMSCF
      common /SCRACH/ STORE(MAXPAR**2)
      dimension XPARAM(MAXPAR), GR(3,NUMATM),
     1DELDIP(3,MAXPAR), TRDIP(3,MAXPAR),LOCOLD(2,MAXPAR)
     2,REDMAS(MAXPAR), SHIFT(6), DIPT(MAXPAR), TRAVEL(MAXPAR)
     3, ROT(3,3), GEOREF(3,MAXPAR)
      character KEYWRD*241, KEYS(241)*1, ELEMNT*2, SIMBOL*10
      logical RESTRT, LINEAR, DEBUG, BARTEL, PRNT, LARGE, LIMSCF
      equivalence (GRAD(1), GR(1,1)), (KEYWRD,KEYS(1))
C
C TEST GEOMETRY TO SEE IF IT IS OPTIMIZED
      TIME2=-1.D9
      call GMETRY(GEO,COORD)
      NVAOLD=NVAR
      do 10 I=1,NVAR
         LOCOLD(1,I)=LOC(1,I)
   10 LOCOLD(2,I)=LOC(2,I)
      NVAR=0
      NDEOLD=NDEP
      NDEP=0
      NUMAT=0
      if (LABELS(1) .NE. 99) NUMAT=1
      do 30 I=2,NATOMS
         if (LABELS(I).EQ.99) go to 30
         if (I.EQ.2) ILIM=1
         if (I.EQ.3) ILIM=2
         if (I.GT.3) ILIM=3
C
C  IS IT A POLYMER?
C
         if (LABELS(I).EQ.107) then
            ILIM=1
         else
            NUMAT=NUMAT+1
         end if
C$DOIT ASIS
         do 20 J=1,ILIM
            NVAR=NVAR+1
            LOC(1,NVAR)=I
            LOC(2,NVAR)=J
   20    XPARAM(NVAR)=GEO(J,I)
   30 continue
C
C   IF A RESTART, THEN TSCF AND TDER WILL BE FAULTY, THEREFORE SET TO -1
C
      TSCF=-1.D0
      TDER=-1.D0
      PRNT=(INDEX(KEYWRD,'RC=') .EQ. 0)
      DEBUG=(INDEX(KEYWRD,'DFORCE') .NE. 0)
      LARGE=(INDEX(KEYWRD,'LARGE') .NE. 0)
      BARTEL=(INDEX(KEYWRD,'NLLSQ') .NE. 0)
      RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0)
      TIME1=SECOND()
      if (RESTRT) then
C
C   CHECK TO SEE if CALCULATION IS IN NLLSQ OR FORCE.
C
         if (BARTEL) go to 50
C
C   CALCULATION IS IN FORCE
C
         go to 80
      end if
*
*      write(6,*) '  !'
*      write(6,*) ' Let the  F O R C E be with you !'  ! J/M dbgng
*
      call COMPFG( XPARAM, .TRUE., ESCF, .TRUE., GRAD, .FALSE.)
      if (PRNT) write(MIOB,'(//10X,''HEAT OF FORMATION ='',F12.6,
     %                            '' KCALS/MOLE'')')   ESCF
      TIME2=SECOND()
      TSCF=TIME2-TIME1
      call COMPFG( XPARAM, .TRUE., ESCF1, .FALSE., GRAD, .TRUE.)
      TIME3=SECOND()
      TDER=TIME3-TIME2
      if (PRNT) write(MIOB,
     %                 '(//10X,''INTERNAL COORDINATE DERIVATIVES'',//3X,
     1 ''NUMBER  ATOM'',2X,''BOND'',9X,''  ANGLE'',10X,''DIHEDRAL'',/)')
      L=0
      IU=0
      do 40 I=1,NATOMS
         if (LABELS(I).EQ.99) go to 40
         L=L+1
         IL=IU+1
         if (I .EQ. 1) IU=IL-1
         if (I .EQ. 2) IU=IL
         if (I .EQ. 3) IU=IL+1
         if (I .GT. 3) IU=IL+2
         if (LABELS(I).EQ.107) IU=IL
         if (PRNT)  write(MIOB,'(I6,4X,A2,F13.6,2F13.6)')
     %                        L,ELEMNT(LABELS(I)),(GRAD(J),J=IL,IU)
   40 continue
C   TEST SUM OF GRADIENTS
      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
      if (PRNT) write(MIOB,'(//10X,''GRADIENT NORM ='',F10.5)') GNORM
      if (GNORM.LT.10.D0) go to 70
      if (INDEX(KEYWRD,' LET ') .NE. 0) then
         write(MIOB,'(///1X,
     %    ''** GRADIENT IS VERY LARGE, BUT SINCE "LET"'',
     %    '' IS USED, CALCULATION WILL continue'')')
         go to 80
      end if
      write(MIOB,'(///1X,''** GRADIENT IS TOO LARGE TO ALLOW '',
     1    ''FORCE MATRIX TO BE CALCULATED, (LIMIT=10) **'',//)')
   50 continue
      do 60 I=1,NVAR
   60 SIMBOL(I)='---'
      write(MIOB,'(//10X,'' GEOMETRY WILL BE OPTIMIZED FIRST'')')
      if (BARTEL) then
         write(MIOB,'(15X,''USING NLLSQ'')')
         call NLLSQ(XPARAM,NVAR)
      else
         write(MIOB,'(15X,''USING FLEPO'')')
         call FLEPO(XPARAM,NVAR,ESCF)
C
C  DID FLEPO USE ALL THE TIME ALLOWED?
C
         if (IFLEPO.EQ.-1) return
      end if
      LIMSCF=.FALSE.
      call COMPFG( XPARAM, .TRUE., ESCF, .TRUE., GRAD, .TRUE.)
      call WRITMO(TIME1,ESCF)
      write(MIOB,'(//10X,''GRADIENT NORM ='',F10.7)') GNORM
      call GMETRY(GEO,COORD)
   70 continue
      do 8 I=1,3
      do 8 J=1,NATOMS
   8  GEOREF(I,J)=GEO(I,J)
C
C NOW TO CALCULATE THE FORCE MATRIX
C
C CHECK OUT SYMMETRY
   80 continue
C
C   NEED TO ENSURE THAT XYZINT WILL WORK CORRECTLY BEFORE call
C   TO DRC.
C
      L=0
      do 90 I=1,NATOMS
         if (LABELS(I).NE.99) then
            L=L+1
            LABELS(L)=LABELS(I)
         end if
   90 continue
      NATOMS=NUMAT
      call XYZINT(COORD,NUMAT,NA,NB,NC,1.D0,GEO)
      call GMETRY(GEO,COORD)
      if (INDEX(KEYWRD,'THERMO').NE.0 .AND.GNORM.GT.1.D0) then
         write(MIOB,'(//30X,''**** WARNING ****'',//
     %    10X,'' GRADIENT IS VERY LARGE FOR A THERMO CALCULATION'',/
     %    10X,'' RESULTS ARE LIKELY TO BE INACCURATE if THERE ARE'')')
         write(MIOB,
     %        '(10X,'' ANY LOW-LYING VIBRATIONS (LESS THAN ABOUT '',
     %        ''400CM-1)'')')
         write(MIOB,'(10X,'' GRADIENT NORM SHOULD BE LESS THAN ABOUT '',
     %          ''0.2 FOR THERMO'',/10X,'' TO GIVE ACCURATE RESULTS'')')
      end if
      if (TSCF.GT.0.D0) then
         write(MIOB,'(//10X,''TIME FOR SCF CALCULATION ='',F8.2)') TSCF
         write(MIOB,'(//10X,''TIME FOR DERIVATIVES     ='',F8.2)') TDER
      end if
      if (NDEP.GT.0) then
         write(MIOB,'(//10X,''SYMMETRY WAS SPECifIED, BUT '',
     %                      ''CANNOT BE USED HERE'')')
         NDEP=0
      end if
      if (PRNT) call AXIS(COORD,NUMAT,A,B,C,WTMOL,2,ROT)
      NVIB=3*NUMAT-6
      if (ABS(C).LT.1.D-20) NVIB=NVIB+1
      if (ID.NE.0)          NVIB=3*NUMAT-3
      if (PRNT) then
         write(MIOB,'(/9X,
     %              ''ORIENTATION OF MOLECULE IN FORCE CALCULATION'')')
         write(MIOB,'(/,4X,''NO.'',7X,''ATOM'',9X,''X'',
     %                  9X,''Y'',9X,''Z'',/)')
      end if
      L=0
      do 100 I=1,NATOMS
         if (LABELS(I) .EQ. 99) go to 100
         L=L+1
         if (PRNT) write(MIOB,'(I6,7X,I3,4X,3F10.4)')
     %                                  L,LABELS(I),(COORD(J,L),J=1,3)
  100 continue
      call FMAT(FMATRX, NVIB, TSCF, TDER, DELDIP,ESCF)
      NA(1)=0
      do 18 I=1,3
      do 18 J=1,NATOMS
   18  GEO(I,J)=GEOREF(I,J)
      if (NVIB.LT.0) then
      NDEP=NDEOLD
      NVAR=0
      return
      end if
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 130 I=1,NUMAT
         IL=IU+1
         IU=IL+2
         IM1=I-1
         JU=0
         do 120 J=1,IM1
            JL=JU+1
            JU=JL+2
            SUM=0.D0
C$DOIT ASIS
            do 110 II=IL,IU
C$DOIT ASIS
               do 110 JJ=JL,JU
  110       SUM=SUM+FMATRX((II*(II-1))/2+JJ)**2
            IJ=IJ+1
  120    STORE(IJ)=SQRT(SUM)
         IJ=IJ+1
  130 STORE(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))
      if (DEBUG) then
         write(MIOB,
     %         '(//10X,'' FULL FORCE MATRIX, INVOKED BY "DFORCE"'')')
         I=-NVAR
         call VECPRT(FMATRX,I)
      end if
      if (PRNT) then
         write(MIOB,'(//10X,'' FORCE MATRIX IN MILLIDYNES/ANGSTROM'')')
         call VECPRT(STORE,NUMAT)
      end if
      L=(NVAR*(NVAR+1))/2
      do 140 I=1,L
  140 STORE(I)=FMATRX(I)
      if (PRNT) call AXIS(COORD,NUMAT,A,B,C,SUM,0,ROT)
      if (PRNT) write(MIOB,'(//10X,''HEAT OF FORMATION ='',F12.6,
     %                            '' KCALS/MOLE'')')  ESCF
      if (LARGE) then
         call FRAME(STORE,NUMAT,0, SHifT)
         call RSP(STORE,NVAR,NVAR,FREQ,CNORML)
         do 150 I=NVIB+1,NVAR
            J=(FREQ(I)+50.D0)*0.01D0
  150    FREQ(I)=FREQ(I)-J*100
         if (PRNT) then
            write(MIOB,'(//10X,''TRIVIAL VIBRATIONS, SHOULD BE ZERO'')')
            write(MIOB,'(/, F9.4,''=TX'',F9.4,''=TY'',F9.4,''=TZ'',
     %             F9.4,''=RX'',F9.4,''=RY'',F9.4,''=RZ'')')
     %                                           (FREQ(I),I=NVIB+1,NVAR)
         write(MIOB,'(//10X,''FORCE CONSTANTS IN MILLIDYNES/ANGSTROM''
     %                        ,'' (= 10**5 DYNES/CM)'',/)')
            write(MIOB,'(8F10.5)') (FREQ(I),I=1,NVIB)
C CONVERT TO WEIGHTED FMAT
            write(MIOB,'(//10X,'' ASSOCIATED EIGENVECTORS'')')
            I=-NVAR
            call MATOUT(CNORML,FREQ,NVIB,I,NVAR)
         end if
      end if
      call FREQCY(FMATRX,FREQ,CNORML,REDMAS,TRAVEL,.TRUE.)
C
C  CALCULATE ZERO POINT ENERGY
C
C
C  THESE CONSTANTS TAKEN FROM HANDBOOK OF CHEMISTRY AND PHYSICS 62ND ED.
C   N AVOGADRO'S NUMBER = 6.022045*10**23
C   H PLANCK'S CONSTANT = 6.626176*10**(-34)JHZ
C   C SPEED OF LIGHT    = 2.99792458*10**10 CM/SEC
C   CONST=0.5*N*H*C/(1000*4.184)
      CONST=1.4295718D-3
      SUM=0.D0
      do 160 I=1,NVAR
  160 SUM=SUM+FREQ(I)
      SUM=SUM*CONST
      if (PRNT)
     1   write(MIOB,'(//10X,'' ZERO POINT ENERGY''
     2              , F12.3,'' KILOCALORIES PER MOLE'')') SUM
      SUMM=0.D0
      do 210 I=1,NVAR
         SUM1=1.D-20
C$DOIT VBEST
         do 170 J=1,NVAR
  170    SUM1=SUM1+CNORML(J+(I-1)*NVAR)**2
         SUM1=1.D0/SQRT(SUM1)
C$DOIT ASIS
         do 180 K=1,3
  180    GRAD(K)=0.D0
C$DOIT ASIS
         do 200 K=1,3
            SUM=0.D0
C$DOIT VBEST
            do 190 J=1,NVAR
  190       SUM=SUM+CNORML(J+(I-1)*NVAR)*DELDIP(K,J)
            SUMM=SUMM+ABS(SUM)
  200    TRDIP(K,I)=SUM*SUM1
         DIPT(I)=SQRT(TRDIP(1,I)**2+TRDIP(2,I)**2+TRDIP(3,I)**2)
  210 continue
      if (PRNT) then
         write(MIOB,'(//3X,'' THE LAST'',I2,'' VIBRATIONS ARE THE'',
     %              '' TRANSLATION AND ROTATION MODES'')') NVAR-NVIB
         write(MIOB,
     %         '(3X,'' THE FIRST THREE OF THESE BEING TRANSLATIONS'',
     %              '' IN X, Y, AND Z, RESPECTIVELY'')')
      end if
      if (PRNT.AND.LARGE) then
         write(MIOB,'(//10X,'' FREQUENCIES, REDUCED MASSES AND '',
     %                      ''VIBRATIONAL DIPOLES''/)')
         NTO6=NVAR/6
         NREM6=NVAR-NTO6*6
         IINC1=-5
         if (NTO6.LT.1) go to 230
         do 220 I=1,NTO6
            write(MIOB,'(/)')
            IINC1=IINC1+6
            IINC2=IINC1+5
            write(MIOB,'(3X,''I'',10I10)') (J,J=IINC1,IINC2)
            write(MIOB,'('' FREQ(I)'',6F10.4,/)')
     %                             (  FREQ(J),J=IINC1,IINC2)
            write(MIOB,'('' MASS(I)'',6F10.5,/)')
     %                             (REDMAS(J),J=IINC1,IINC2)
            write(MIOB,'('' DIPX(I)'',6F10.5)')
     %                            (TRDIP(1,J),J=IINC1,IINC2)
            write(MIOB,'('' DIPY(I)'',6F10.5)')
     %                            (TRDIP(2,J),J=IINC1,IINC2)
            write(MIOB,'('' DIPZ(I)'',6F10.5,/)')
     %                            (TRDIP(3,J),J=IINC1,IINC2)
            write(MIOB,'('' DIPT(I)'',6F10.5)')
     %                               (DIPT(J),J=IINC1,IINC2)
  220    continue
  230    continue
         if (NREM6.LT.1) go to 240
         write(MIOB,'(/)')
         IINC1=IINC1+6
         IINC2=IINC1+(NREM6-1)
         write(MIOB,'(3X,''I'',10I10)')     (J,J=IINC1,IINC2)
         write(MIOB,'('' FREQ(I)'',6F10.4)')
     %                                (FREQ(J),J=IINC1,IINC2)
         write(MIOB,'(/,'' MASS(I)'',6F10.5)')
     %                              (REDMAS(J),J=IINC1,IINC2)
         write(MIOB,'(/,'' DIPX(I)'',6F10.5)')
     %                             (TRDIP(1,J),J=IINC1,IINC2)
         write(MIOB,'('' DIPY(I)'',6F10.5)')
     %                             (TRDIP(2,J),J=IINC1,IINC2)
         write(MIOB,'('' DIPZ(I)'',6F10.5)')
     %                             (TRDIP(3,J),J=IINC1,IINC2)
         write(MIOB,'(/,'' DIPT(I)'',6F10.5)')
     %                                (DIPT(J),J=IINC1,IINC2)
  240    continue
      end if
      if (PRNT) then
         write(MIOB,'(//10X,'' NORMAL COORDINATE ANALYSIS'')')
         I=-NVAR
         call MATOUT(CNORML,FREQ,NVAR,I,NVAR)
      end if
C
C   CARRY OUT IRC if REQUESTED.
C
      if (INDEX(KEYWRD,'IRC')+INDEX(KEYWRD,'DRC').NE.0) then
         do 250 I=1,NVAR
            LOC(1,I)=0
  250    LOC(2,I)=0
         NVAR=NVAOLD
         do 260 I=1,NVAR
            LOC(1,I)=LOCOLD(1,I)
  260    LOC(2,I)=LOCOLD(2,I)
         call XYZINT(COORD,NUMAT,NA,NB,NC,1.D0,GEO)
         LAST=1
         call DRC(CNORML,FREQ)
      NA(1)=0
      NDEP=NDEOLD
      NVAR=0
      do 28 I=1,3
      do 28 J=1,NATOMS
   28  GEO(I,J)=GEOREF(I,J)
         return
      end if
      call FREQCY(FMATRX,FREQ,CNORML,DELDIP,DELDIP,.FALSE.)
      write(MIOB,'(//10X,'' MASS-WEIGHTED COORDINATE ANALYSIS'')')
      I=-NVAR
      call MATOUT(CNORML,FREQ,NVAR,I,NVAR)
      call ANAVIB(COORD,FREQ,DIPT,NVAR,CNORML,STORE,
     1FMATRX,TRAVEL,REDMAS)
      if (INDEX(KEYWRD,'THERMO').NE.0) then
         call GMETRY(GEO,COORD)
         I=INDEX(KEYWRD,' ROT')
         if (I.NE.0) then
            SYM = READA(KEYWRD,I)
         else
            SYM=1
         end if
         LINEAR=(ABS(A*B*C) .LT. 1.D-10)
         I=INDEX(KEYWRD,' TRANS')
C
C   "I" IS GOING TO MARK THE BEGINNING OF THE GENUINE VIBRATIONS.
C
         if (I.NE.0) then
            I=INDEX(KEYWRD,' TRANS=')
            if (I.NE.0) then
               I=1+READA(KEYWRD,I)
               J=NVIB-I+1
               write(MIOB,
     %               '(//1X,''THE LOWEST'',I3,'' VIBRATIONS ARE NOT'',
     %               /,'' TO BE USED IN THE THERMO CALCULATION'')') I-1
            else
               write(MIOB,'(//10X,''SYSTEM IS A TRANSITION STATE'')')
               I=2
               J=NVIB-1
            end if
         else
            write(MIOB,'(//10X,''SYSTEM IS A GROUND STATE'')')
            I=1
            J=NVIB
         end if
         call THERMO(A,B,C,LINEAR,SYM,WTMOL,FREQ(I),J,ESCF)
      end if
      NA(1)=0
      NVAR=0
      NDEP=NDEOLD
      do 38 I=1,3
      do 38 J=1,NATOMS
   38  GEO(I,J)=GEOREF(I,J)
      return
      end
