      SUBROUTINE FORCE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'mopac.inc'
      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 /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 /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)
      CHARACTER KEYWRD*80, KEYS(80)*1, ELEMNT*2
      LOGICAL RESTRT, LINEAR, DEBUG, BARTEL, PRNT, LARGE
      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
      NUMAT=0
      IF(LABELS(1) .NE. 99) NUMAT=1
      DO 30 I=2,NATOMS
         IF(LABELS(I).EQ.99.OR.LABELS(I).EQ.107) GOTO 30
         NUMAT=NUMAT+1
         IF(I.EQ.2)ILIM=1
         IF(I.EQ.3)ILIM=2
         IF(I.GT.3)ILIM=3
         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=ZECOND()
      IF (RESTRT) THEN
C
C   CHECK TO SEE IF CALCULATION IS IN NLLSQ OR FORCE.
C
         IF(BARTEL)GOTO 50
C
C   CALCULATION IS IN FORCE
C
         GOTO 70
      ENDIF
      CALL COMPFG( XPARAM, .TRUE., ESCF, .TRUE., GRAD, .FALSE.)
      IF(PRNT)WRITE(6,'(//10X,''HEAT OF FORMATION ='',F12.6,
     1'' KCALS/MOLE'')')ESCF
      TIME2=ZECOND()
      TSCF=TIME2-TIME1
      CALL COMPFG( XPARAM, .TRUE., ESCF1, .FALSE., GRAD, .TRUE.)
      TIME3=ZECOND()
      TDER=TIME3-TIME2
      IF(PRNT)WRITE(6,'(//10X,''INTERNAL COORDINATE DERIVATIVES'',//3X,
     1''ATOM  AT. NO.'',2X,''BOND'',9X,''ANGLE'',8X,''DIHEDRAL'',/)')
      L=0
      IU=0
      DO 40 I=1,NATOMS
         IF(LABELS(I).EQ.99) GOTO 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(PRNT)WRITE(6,'(I6,4X,A2,F13.6,2F10.6)')
     1L,ELEMNT(LABELS(I)),(GRAD(J),J=IL,IU)
   40 CONTINUE
C   TEST SUM OF GRADIENTS
      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
      IF(PRNT)WRITE(6,'(//10X,''GRADIENT NORM ='',F10.5)') GNORM
      IF(GNORM.LT.10.D0) GOTO 60
      IF(INDEX(KEYWRD,' LET ') .NE. 0) THEN
         WRITE(6,'(///1X,''** GRADIENT IS VERY LARGE, BUT SINCE "LET"'',
     1'' IS USED, CALCULATION WILL CONTINUE'')')
         GOTO 70
      ENDIF
      WRITE(6,'(///1X,''** GRADIENT IS TOO LARGE TO ALLOW '',
     1    ''FORCE MATRIX TO BE CALCULATED, (LIMIT=10) **'',//)')
   50 CONTINUE
      WRITE(6,'(//10X,'' GEOMETRY WILL BE OPTIMIZED FIRST'')')
      IF(BARTEL) THEN
         WRITE(6,'(15X,''USING NLLSQ'')')
         CALL NLLSQ(XPARAM,NVAR)
      ELSE
         WRITE(6,'(15X,''USING FLEPO'')')
         CALL FLEPO(XPARAM,NVAR,ESCF)
      ENDIF
      CALL COMPFG( XPARAM, .TRUE., ESCF, .TRUE., GRAD, .TRUE.)
      CALL VWRITE(TIME1,ESCF)
      WRITE(6,'(//10X,''GRADIENT NORM ='',F10.7)') GNORM
      CALL GMETRY(GEO,COORD)
   60 CONTINUE
C
C NOW TO CALCULATE THE FORCE MATRIX
C
C CHECK OUT SYMMETRY
   70 CONTINUE
C
C   NEED TO ENSURE THAT XYZINT WILL WORK CORRECTLY BEFORE CALL
C   TO DRC.
C
      L=0
      DO 80 I=1,NATOMS
         IF(LABELS(I).NE.99)THEN
            L=L+1
            LABELS(L)=LABELS(I)
         ENDIF
   80 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(6,'(//30X,''**** WARNING ****'',//
     110X,'' GRADIENT IS VERY LARGE FOR A THERMO CALCULATION'',/
     210X,'' RESULTS ARE LIKELY TO BE INACCURATE IF THERE ARE'')')
         WRITE(6,'(10X,'' ANY LOW-LYING VIBRATIONS (LESS THAN ABOUT ''
     1,''400CM-1)'')')
         WRITE(6,'(10X,'' GRADIENT NORM SHOULD BE LESS THAN ABOUT '',
     1''0.2 FOR THERMO'',/10X,'' TO GIVE ACCURATE RESULTS'')')
      ENDIF
      IF(TSCF.GT.0.D0) THEN
         WRITE(6,'(//10X,''TIME FOR SCF CALCULATION ='',F8.2)')TSCF
         WRITE(6,'(//10X,''TIME FOR DERIVATIVES     ='',F8.2)')TDER
      ENDIF
      IF(NDEP.GT.0) THEN
         WRITE(6,'(//10X,''SYMMETRY WAS SPECIFIED, BUT '',
     1''CANNOT BE USED HERE'')')
         NDEP=0
      ENDIF
      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(6,'(/9X,''ORIENTATION OF MOLECULE IN FORCE CALCULATION'')
     1')
         WRITE(6,'(/,4X,''NO.'',7X,''ATOM'',9X,''X'',
     19X,''Y'',9X,''Z'',/)')
      ENDIF
      L=0
      DO 90 I=1,NATOMS
         IF(LABELS(I) .EQ. 99) GOTO 90
         L=L+1
         IF(PRNT)WRITE(6,'(I6,7X,I3,4X,3F10.4)')
     1    L,LABELS(I),(COORD(J,L),J=1,3)
   90 CONTINUE
      CALL FMAT(FMATRX, NVIB, TSCF, TDER, DELDIP,ESCF)
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 120 I=1,NUMAT
         IL=IU+1
         IU=IL+2
         IM1=I-1
         JU=0
         DO 110 J=1,IM1
            JL=JU+1
            JU=JL+2
            SUM=0.D0
            DO 100 II=IL,IU
               DO 100 JJ=JL,JU
  100       SUM=SUM+FMATRX((II*(II-1))/2+JJ)**2
            IJ=IJ+1
  110    STORE(IJ)=SQRT(SUM)
         IJ=IJ+1
  120 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(6,'(//10X,'' FULL FORCE MATRIX, INVOKED BY "DFORCE"'')')
         I=-NVAR
         CALL VECPRT(FMATRX,I)
      ENDIF
      IF(PRNT)THEN
         WRITE(6,'(//10X,'' FORCE MATRIX IN MILLIDYNES/ANGSTROM'')')
         CALL VECPRT(STORE,NUMAT)
      ENDIF
      L=(NVAR*(NVAR+1))/2
      DO 130 I=1,L
  130 STORE(I)=FMATRX(I)
      IF(PRNT) CALL AXIS(COORD,NUMAT,A,B,C,SUM,0,ROT)
      IF(PRNT)WRITE(6,'(//10X,''HEAT OF FORMATION ='',F12.6,
     1'' KCALS/MOLE'')')ESCF
      IF(LARGE)THEN
         CALL FRAME(STORE,NUMAT,0, SHIFT)
         CALL RSP(STORE,NVAR,NVAR,FREQ,CNORML)
         DO 140 I=NVIB+1,NVAR
            J=(FREQ(I)+50.D0)*0.01D0
  140    FREQ(I)=FREQ(I)-J*100
         IF(PRNT)THEN
            WRITE(6,'(//10X,''TRIVIAL VIBRATIONS, SHOULD BE ZERO'')')
            WRITE(6,'(/, F9.4,''=TX'',F9.4,''=TY'',F9.4,''=TZ'',
     1             F9.4,''=RX'',F9.4,''=RY'',F9.4,''=RZ'')')
     2(FREQ(I),I=NVIB+1,NVAR)
            WRITE(6,'(//10X,''FORCE CONSTANTS IN MILLIDYNES/ANGSTROM''
     1,'' (= 10**5 DYNES/CM)'',/)')
            WRITE(6,'(8F10.5)')(FREQ(I),I=1,NVIB)
C CONVERT TO WEIGHTED FMAT
            WRITE(6,'(//10X,'' ASSOCIATED EIGENVECTORS'')')
            I=-NVAR
            CALL MATOUT(CNORML,FREQ,NVIB,I,NVAR)
         ENDIF
      ENDIF
      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 150 I=1,NVAR
  150 SUM=SUM+FREQ(I)
      SUM=SUM*CONST
      IF(PRNT)
     1WRITE(6,'(//10X,'' ZERO POINT ENERGY''
     2, F12.3,'' KILOCALORIES PER MOLE'')')SUM
      SUMM=0.D0
      DO 200 I=1,NVAR
         SUM1=1.D-20
         DO 160 J=1,NVAR
  160    SUM1=SUM1+CNORML(J+(I-1)*NVAR)**2
         SUM1=1.D0/SQRT(SUM1)
         DO 170 K=1,3
  170    GRAD(K)=0.D0
         DO 190 K=1,3
            SUM=0.D0
            DO 180 J=1,NVAR
  180       SUM=SUM+CNORML(J+(I-1)*NVAR)*DELDIP(K,J)
            SUMM=SUMM+ABS(SUM)
  190    TRDIP(K,I)=SUM*SUM1
         DIPT(I)=SQRT(TRDIP(1,I)**2+TRDIP(2,I)**2+TRDIP(3,I)**2)
  200 CONTINUE
      IF(PRNT)THEN
         WRITE(6,'(//3X,'' THE LAST'',I2,'' VIBRATIONS ARE THE'',
     1'' TRANSLATION AND ROTATION MODES'')')NVAR-NVIB
         WRITE(6,'(3X,'' THE FIRST THREE OF THESE BEING TRANSLATIONS'',
     1'' IN X, Y, AND Z, RESPECTIVELY'')')
      ENDIF
      IF(PRNT.AND.LARGE)THEN
         WRITE(6,'(//10X,'' FREQUENCIES, REDUCED MASSES AND '',
     1''VIBRATIONAL DIPOLES''/)')
         NTO6=NVAR/6
         NREM6=NVAR-NTO6*6
         IINC1=-5
         IF (NTO6.LT.1) GO TO 220
         DO 210 I=1,NTO6
            WRITE (6,'(/)')
            IINC1=IINC1+6
            IINC2=IINC1+5
            WRITE (6,'(3X,''I'',10I10)') (J,J=IINC1,IINC2)
            WRITE (6,'('' FREQ(I)'',6F10.4,/)') (FREQ(J),J=IINC1,IINC2)
            WRITE (6,'('' MASS(I)'',6F10.5,/)') (REDMAS(J),J=IINC1,IINC2
     1)
            WRITE (6,'('' DIPX(I)'',6F10.5)') (TRDIP(1,J),J=IINC1,IINC2)
            WRITE (6,'('' DIPY(I)'',6F10.5)') (TRDIP(2,J),J=IINC1,IINC2)
            WRITE (6,'('' DIPZ(I)'',6F10.5,/)') (TRDIP(3,J),J=IINC1,IINC
     12)
            WRITE (6,'('' DIPT(I)'',6F10.5)')
     1   (DIPT(J),J=IINC1,IINC2)
  210    CONTINUE
  220    CONTINUE
         IF (NREM6.LT.1) GO TO 230
         WRITE (6,'(/)')
         IINC1=IINC1+6
         IINC2=IINC1+(NREM6-1)
         WRITE (6,'(3X,''I'',10I10)') (J,J=IINC1,IINC2)
         WRITE (6,'('' FREQ(I)'',6F10.4)') (FREQ(J),J=IINC1,IINC2)
         WRITE (6,'(/,'' MASS(I)'',6F10.5)') (REDMAS(J),J=IINC1,IINC2)
         WRITE (6,'(/,'' DIPX(I)'',6F10.5)') (TRDIP(1,J),J=IINC1,IINC2)
         WRITE (6,'('' DIPY(I)'',6F10.5)') (TRDIP(2,J),J=IINC1,IINC2)
         WRITE (6,'('' DIPZ(I)'',6F10.5)') (TRDIP(3,J),J=IINC1,IINC2)
         WRITE (6,'(/,'' DIPT(I)'',6F10.5)')
     1   (DIPT(J),J=IINC1,IINC2)
  230    CONTINUE
      ENDIF
      WRITE(6,'(//10X,'' NORMAL COORDINATE ANALYSIS'')')
      I=-NVAR
      CALL MATOUT(CNORML,FREQ,NVAR,I,NVAR)
C
C   CARRY OUT IRC IF REQUESTED.
C
      IF(INDEX(KEYWRD,'IRC')+INDEX(KEYWRD,'DRC').NE.0)THEN
         DO 240 I=1,NVAR
            LOC(1,I)=0
  240    LOC(2,I)=0
         NVAR=NVAOLD
         DO 250 I=1,NVAR
            LOC(1,I)=LOCOLD(1,I)
  250    LOC(2,I)=LOCOLD(2,I)
         CALL XYZINT(COORD,NUMAT,NA,NB,NC,1.D0,GEO)
         LAST=1
         CALL DRC(CNORML,FREQ)
         CALL WRAPUP()
      ENDIF
      CALL FREQCY(FMATRX,FREQ,CNORML,DELDIP,DELDIP,.FALSE.)
      WRITE(6,'(//10X,'' MASS-WEIGHTED COORDINATE ANALYSIS'')')
      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
         ENDIF
         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(6,'(//1X,''THE LOWEST'',I3,'' VIBRATIONS ARE NOT'',
     1/,'' TO BE USED IN THE THERMO CALCULATION'')')I-1
            ELSE
               WRITE(6,'(//10X,''SYSTEM IS A TRANSITION STATE'')')
               I=2
               J=NVIB-1
            ENDIF
         ELSE
            WRITE(6,'(//10X,''SYSTEM IS A GROUND STATE'')')
            I=1
            J=NVIB
         ENDIF
         CALL THERMO(A,B,C,LINEAR,SYM,WTMOL,FREQ(I),J,ESCF)
      ENDIF
      RETURN
      END
