      SUBROUTINE FMAT(FMATRX, TSCF, TDER, DELDIP, HEAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION FMATRX(*), DELDIP(3,*)
***********************************************************************
*
*  VALUE CALCULATES THE SECOND-ORDER OF THE ENERGY WITH
*        RESPECT TO THE CARTESIAN COORDINATES I AND J AND PLACES IT
*        IN FMATRX
*
*  ON INPUT NATOMS  = NUMBER OF ATOMS IN THE SYSTEM.
*           XPARAM  = INTERNAL COORDINATES OF MOLECULE STORED LINEARLY
*
*  VARIABLES USED
*           COORDL  = ARRAY OF CARTESIAN COORDINATES, STORED LINEARLY.
*           I       = INDEX OF CARTESIAN COORDINATE.
*           J       = INDEX OF CARTESIAN COORDINATE.
*
*  ON OUTPUT FMATRX = SECOND DERIVATIVE OF THE ENERGY WITH RESPECT TO
*                    CARTESIAN COORDINATES I AND J.
***********************************************************************
      COMMON /KEYWRD/ KEYWRD
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY, DUMY(MAXPAR)
      COMMON /DENSTY/ P(MPACK),PDUMY(2,MPACK)
      COMMON /TIME  / TIME0
      COMMON /CORE  / CORE(107)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /COORD / COORD
      COMMON /SCRACH/ EVECS(MAXPAR*MAXPAR)
      DIMENSION COLD(MAXPAR), GRAD(MAXPAR), COORD(3,NUMATM),
     1GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3)
     2, EIGS(MAXPAR),
     3 FCONST(MAXPAR)
      CHARACTER*80 KEYWRD
      LOGICAL DEBUG, DERIV, RESTRT
      CHARACTER SPACE*1, DOT*1, ZERO*1, NINE*1, CH*1
      EQUIVALENCE (COORD(1,1),COORDL(1))
      DATA SPACE,DOT,ZERO,NINE /' ','.','0','9'/
      DATA FACT/6.95125D-3/
C
C    FACT IS THE CONVERSION FACTOR FROM KCAL/MOLE TO ERGS
C
C SET UP CONSTANTS AND FLAGS
      NA(1)=99
C
C  SET UP THE VARIABLES IN XPARAM ANDLOC,THESE ARE IN CARTESIAN COORDINA
C
      NUMAT=0
      DO 10 I=1,NATOMS
         IF(LABELS(I).NE.99.AND.LABELS(I).NE.107) THEN
            NUMAT=NUMAT+1
            LABELS(NUMAT)=LABELS(I)
         ENDIF
   10 CONTINUE
      NATOMS=NUMAT
C
C   THIS IS A QUICK, IF CLUMSY, WAY TO CALCULATE NUMAT, AND TO REMOVE TH
C   DUMMY ATOMS FROM THE ARRAY LABELS.
C
      NVAR=0
      DO 20 I=1,NUMAT
         DO 20 J=1,3
            NVAR=NVAR+1
            LOC(1,NVAR)=I
            LOC(2,NVAR)=J
   20 CONTINUE
      LIN=(NVAR*(NVAR+1))/2
      DO 30 I=1,LIN
   30 FMATRX(I)=0.D0
      RESTRT =(INDEX(KEYWRD,'RESTART') .NE. 0)
      IF(INDEX(KEYWRD,'NLLSQ') .NE. 0) RESTRT=.FALSE.
      DEBUG =(INDEX(KEYWRD,'FMAT') .NE. 0)
      DERIV=(NCLOSE .EQ. NOPEN .AND. INDEX(KEYWRD,'C.I.') .EQ. 0)
      WRITE(6,'(//4X,''FIRST DERIVATIVES WILL BE USED IN THE''
     1,'' CALCULATION OF SECOND DERIVATIVES'')')
      TIME=3600
      I=INDEX(KEYWRD,' T=')
      IF(I.NE.0) THEN
         TIM=READA(KEYWRD,I)
         DO 40 J=I+3,80
            CH=KEYWRD(J:J)
            IF( CH .NE. DOT .AND. (CH .LT. ZERO .OR. CH .GT. NINE)) THEN
               IF( CH .EQ. 'M') TIM=TIM*60
               GOTO 50
            ENDIF
   40    CONTINUE
   50    TIME=TIM
         WRITE(6,'(/10X,''TIME DEFINED FOR THIS STEP ='',F19.2,
     1    '' SECONDS'')')TIME
      ELSE
         WRITE(6,'(/10X,''DEFAULT TIME OF'',F8.2,
     1    '' SECONDS ALLOCATED FOR THIS STEP'')')TIME
      ENDIF
      TLEFT  = TIME-SECOND()+TIME0
      IF(RESTRT) THEN
         DO 60 I=1,NVAR
   60    COLD(I)=COORDL(I)
         ISTART = 0
         I=0
         CALL FORSAV(TOTIME,DELDIP,ISTART,I,FMATRX, COORD, NVAR,HEAT,
     1                EVECS,JSTART,FCONST)
         KOUNTF=(ISTART*(ISTART+1))/2
         ISTART=ISTART+1
         JSTART=JSTART+1
         TIME2 = SECOND()
         IF(ISTART.GT.NVAR) GOTO 160
      ELSE
         TOTIME=0.D0
         IF (TSCF.GT.0.D0)TLEFT=TLEFT-TSCF-TDER
         ISTART=1
      ENDIF
      CONST =60
      DELTA =0.5D0/CONST
C CALCULATE FMATRX
      IF(ISTART.GT.1) THEN
         ESTIME=(NVAR-ISTART+1)*TOTIME/(ISTART-1.D0)
      ELSE
         ESTIME=NVAR*(TSCF+TDER)*2.D0
      ENDIF
      IF(TSCF.GT.0)
     1WRITE(6,'(/10X,''ESTIMATED TIME TO COMPLETE CALCULATION =''
     2,F9.2,'' SECONDS'')')ESTIME
      IF(RESTRT) THEN
         IF(ISTART.LE.NVAR)
     1    WRITE(6,'(/10X,''STARTING AGAIN AT LINE'',18X,I4)')ISTART
         WRITE(6,'(/10X,''TIME USED UP TO RESTART ='',F22.2)')TOTIME
      ENDIF
      LU=KOUNTF
      TIME1 = SECOND()
      NUMAT=NVAR/3
      DO 120 I=ISTART,NVAR
         TIME2 = SECOND()
         COORDL(I)=COORDL(I)+DELTA*0.5
         GROLD(1)=100.D0
         CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.)
         CALL CHRGE(P,Q)
         DO 70 II=1,NUMAT
   70    Q(II)=CORE(LABELS(II))-Q(II)
         SUM = DIPOLE(P,Q,COORDL,DELDIP(1,I))
         COORDL(I)=COORDL(I)-DELTA
         GRAD(1)=100.D0
         CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .TRUE.)
         CALL CHRGE(P,Q)
         DO 80 II=1,NUMAT
   80    Q(II)=CORE(LABELS(II))-Q(II)
         SUM = DIPOLE(P,Q,COORDL,DEL2)
         DO 90 II=1,3
   90    DELDIP(II,I)=(DELDIP(II,I)-DEL2(II))*CONST
         LL=LU+1
         LU=LL+I-1
         L=0
         DO 100 KOUNTF=LL,LU
            L=L+1
            FMATRX(KOUNTF)=FMATRX(KOUNTF)+
     1         ((GROLD(L)-GRAD(L)))
     2          *CONST*FACT*0.5D0
  100    CONTINUE
         L=L-1
         DO 110 K=I,NVAR
            L=L+1
            KK=(K*(K-1))/2+I
            FMATRX(KK)=FMATRX(KK)+
     1         ((GROLD(L)-GRAD(L)))
     2          *CONST*FACT*0.5D0
  110    CONTINUE
         COORDL(I)=COORDL(I)+DELTA*0.5D0
         TIME3 = SECOND()
         TSTEP=TIME3-TIME2
         TOTIME= TOTIME+TSTEP
         TLEFT= TLEFT-TSTEP
         WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL ='',
     1F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
         IF(DERIV) THEN
            ESTIM = TOTIME/I
         ELSE
            ESTIM = TOTIME*2.D0/I
         ENDIF
         IF(I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN
            WRITE(6,'(//10X,''- - - - -  TIME  LIMIT - - - - -'')')
            WRITE(6,'(/10X,'' POINT REACHED ='',I4)')I
            WRITE(6,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')')
            WRITE(6,'(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2,
     1'' SECONDS'')')ESTIM
            JSTART=1
            II=I
            CALL FORSAV(TOTIME,DELDIP,II,NVAR,FMATRX, COORD,NVAR,HEAT,
     1                EVECS,JSTART,FCONST)
         ENDIF
  120 CONTINUE
C#      CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT,
C#     +                EVECS,JSTART,FCONST)
      IF(DERIV) GOTO 250
      WRITE(6,'(//10X,'' STARTING TO CALCULATE FORCE CONSTANTS'',/)')
      CALL FRAME(FMATRX,NUMAT,0,EIGS)
      CALL HQRII(FMATRX,NVAR,NVAR,EIGS,EVECS)
      IF(DEBUG) THEN
         WRITE(6,'(''   EIGENVECTORS FROM FIRST CALCULATION'')')
         CALL MATOUT(EVECS,EIGS,NVAR,NVAR,NVAR)
      ENDIF
      L=0
      NREAL=NVAR-6
      DO 140 I=1,NVAR
         DO 140 J=1,I
            L=L+1
            SUM=0.D0
            DO 130 K=1,NREAL
               K1=(K-1)*NVAR+I
               K2=(K-1)*NVAR+J
  130       SUM=SUM+EVECS(K1)*EIGS(K)*EVECS(K2)
  140 FMATRX(L)=SUM
      CALL FRAME(FMATRX,NUMAT,0,EIGS)
      CALL HQRII(FMATRX,NVAR,NVAR,EIGS,EVECS)
      JSTART=1
      DO 150 I=1,NVAR
  150 COLD(I)=COORDL(I)
  160 IF(DERIV) GOTO 250
      DELTA=0.025D0
      L=(JSTART-1)*NVAR
      DO 210 ILOOP=JSTART,NVAR
         J=L
         DO 170 I=1,NVAR
            J=J+1
  170    COORDL(I)=COLD(I)+EVECS(J)*DELTA
         CALL COMPFG(COORDL, .TRUE., HEATA, .TRUE., GRAD, .FALSE.)
         HEATA=HEATA-HEAT
         J=L
         DO 180 I=1,NVAR
            J=J+1
  180    COORDL(I)=COLD(I)-EVECS(J)*DELTA
         CALL COMPFG(COORDL, .TRUE., HEATB, .TRUE., GRAD, .FALSE.)
         HEATB=HEATB-HEAT
         J=L
         DO 190 I=1,NVAR
            J=J+1
  190    COORDL(I)=COLD(I)+EVECS(J)*DELTA*2
         CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .FALSE.)
         HEATAA=HEATAA-HEAT
         J=L
         DO 200 I=1,NVAR
            J=J+1
  200    COORDL(I)=COLD(I)-EVECS(J)*DELTA*2
         CALL COMPFG(COORDL, .TRUE., HEATBB, .TRUE., GRAD, .FALSE.)
         HEATBB=HEATBB-HEAT
         SUM=( (HEATA+HEATB)*16 - (HEATAA+HEATBB) )/12.D0
     1/DELTA*FACT/DELTA*0.5D0
         FCONST(ILOOP)=SUM*0.5D0
         L=L+NVAR
         TIME3 = SECOND()
         TSTEP=TIME3-TIME2
         TIME2=TIME3
         TOTIME= TOTIME+TSTEP
         TLEFT= TLEFT-TSTEP
         WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL ='',
     1F10.2,'' TIME LEFT:'',F10.2)')ILOOP,TSTEP,TOTIME,TLEFT
         ESTIM = TSTEP*5.D0
C
C    5.0 IS A SAFETY FACTOR
C
         IF(ILOOP.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN
            WRITE(6,'(//10X,''- - - - -  TIME  LIMIT - - - - -'')')
            WRITE(6,'(/10X,'' POINT REACHED ='',I4)')ILOOP
            WRITE(6,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')')
            WRITE(6,'(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2,
     1'' SECONDS'')')ESTIM
            IFOR=ILOOP
            IX=NVAR+2
*
* VALUE OF IX IS NOT IMPORTANT. SHOULD NOT BE 0 OR NVAR
*
            CALL FORSAV(TOTIME,DELDIP,IX,NVAR,FMATRX, COORD,NVAR,HEAT,
     1                EVECS,IFOR,FCONST)
         ENDIF
  210 CONTINUE
      L=0
      DO 230 I=1,NVAR
         DO 230 J=1,I
            L=L+1
            SUM=0.D0
            DO 220 K=1,NVAR
               K1=(K-1)*NVAR+I
               K2=(K-1)*NVAR+J
  220       SUM=SUM+EVECS(K1)*FCONST(K)*EVECS(K2)
  230 FMATRX(L)=SUM*2.D0
      DO 240 I=1,NVAR
  240 COORDL(I)=COLD(I)
  250 CONTINUE
      IF(ISTART.LE.NVAR .AND. INDEX(KEYWRD,'ISOTOPE') .NE. 0)
     1CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT,
     2                EVECS,ILOOP,FCONST)
      RETURN
      END
