      SUBROUTINE FMAT(FMATRX, NREAL, 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 /TIMDMP/ TLEFT, TDUMP
      COMMON /ATMASS/ ATMASS(NUMATM)
      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(3,NUMATM)
      COMMON /NLLCOM/ EVECS(MAXPAR*MAXPAR),BMAT(MAXPAR,MAXPAR*2)
      DIMENSION COLD(MAXPAR), GRAD(MAXPAR),
     1GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3), G2OLD(MAXPAR)
     2, EIGS(MAXPAR), G2RAD(MAXPAR),
     3 FCONST(MAXPAR)
      CHARACTER*241 KEYWRD
      SAVE  SPACE, DOTT, ZERO, NINE, FACT
      LOGICAL DEBUG, RESTRT, PRNT, RESFIL, ANALYT, PRECIS, BIG, LOG
      CHARACTER SPACE*1, DOTT*1, ZERO*1, NINE*1
      EQUIVALENCE (COORD(1,1),COORDL(1))
      DATA SPACE,DOTT,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
C$DOIT ASIS
      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
C   THE DUMMY ATOMS FROM THE ARRAY LABELS.
C
      NVAR=NUMAT*3
      DO 20 I=1,NUMAT
         LOC(1,(I-1)*3+1)=I
         LOC(2,(I-1)*3+1)=1
C
         LOC(1,(I-1)*3+2)=I
         LOC(2,(I-1)*3+2)=2
C
         LOC(1,(I-1)*3+3)=I
         LOC(2,(I-1)*3+3)=3
   20 CONTINUE
      LIN=(NVAR*(NVAR+1))/2
      DO 30 I=1,LIN
   30 FMATRX(I)=0.D0
      PRNT   =(INDEX(KEYWRD,'IRC=') .EQ. 0)
      LOG    =(INDEX(KEYWRD,'NOLOG') .EQ. 0)
      PRECIS =(INDEX(KEYWRD,'PREC') .NE. 0)
      ANALYT =(INDEX(KEYWRD,'ANALYT') .NE. 0)
      RESTRT =(INDEX(KEYWRD,'RESTART') .NE. 0)
      IF(INDEX(KEYWRD,'NLLSQ') .NE. 0) RESTRT=.FALSE.
      DEBUG =(INDEX(KEYWRD,'FMAT') .NE. 0)
      BIG    =(INDEX(KEYWRD,'LARGE') .NE. 0 .AND. DEBUG)
      IF(PRNT)WRITE(6,'(//4X,''FIRST DERIVATIVES WILL BE USED IN THE''
     1,'' CALCULATION OF SECOND DERIVATIVES'')')
      TLAST=TLEFT
      RESFIL=.FALSE.
      IF(RESTRT) THEN
         DO 40 I=1,NVAR
   40    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()
      ELSE
         KOUNTF=0
         TOTIME=0.D0
         IF (TSCF.GT.0.D0)TLEFT=TLEFT-TSCF-TDER
         ISTART=1
      ENDIF
C CALCULATE FMATRX
      IF(ISTART.GT.1) THEN
         ESTIME=(NVAR-ISTART+1)*TOTIME/(ISTART-1.D0)
      ELSE
         ESTIME=NVAR*(TSCF+TDER)*2.D0
         IF (PRECIS) ESTIME=ESTIME*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 50 I=1,NVAR
   50 EIGS(I)=0.D0
      DO 120 I=ISTART,NVAR
         TIME2 = SECOND()
         DELTA=1.D0/120.D0
         IF(PRECIS)THEN
C
C   DETERMINE A GOOD STEP SIZE
C
            G2OLD(1)=100.D0
            COORDL(I)=COORDL(I)+DELTA
            CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
            COORDL(I)=COORDL(I)-DELTA
            DELTA=DELTA*10.D0/SQRT(DOT(G2OLD,G2OLD,NVAR))
C
C   CONSTRAIN DELTA TO A 'REASONABLE' VALUE
C
            DELTA=MIN(0.05D0,MAX(0.005D0,DELTA))
            IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
     1A
            G2OLD(1)=100.D0
            COORDL(I)=COORDL(I)+DELTA
            CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
            IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +1.0*DELTA',
     1SQRT(DOT(G2OLD,G2OLD,NVAR))
            COORDL(I)=COORDL(I)-DELTA*2.D0
            G2RAD(1)=100.D0
            CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., G2RAD, .TRUE.)
            COORDL(I)=COORDL(I)+DELTA
            IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -1.0*DELTA',
     1SQRT(DOT(G2RAD,G2RAD,NVAR))
         ELSE
            IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
     1A
         ENDIF
         COORDL(I)=COORDL(I)+0.5D0*DELTA
         GROLD(1)=100.D0
         CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.)
         IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +0.5*DELTA',
     1SQRT(DOT(GROLD,GROLD,NVAR))
         CALL CHRGE(P,Q)
         DO 60 II=1,NUMAT
   60    Q(II)=CORE(LABELS(II))-Q(II)
         SUM = DIPOLE(P,Q,COORDL,DELDIP(1,I),0)
         COORDL(I)=COORDL(I)-DELTA
         GRAD(1)=100.D0
         CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .TRUE.)
         IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -0.5*DELTA',
     1SQRT(DOT(GRAD,GRAD,NVAR))
         CALL CHRGE(P,Q)
         DO 70 II=1,NUMAT
   70    Q(II)=CORE(LABELS(II))-Q(II)
         SUM = DIPOLE(P,Q,COORDL,DEL2,0)
         COORDL(I)=COORDL(I)+DELTA*0.5D0
         DELDIP(1,I)=(DELDIP(1,I)-DEL2(1))*0.5D0/DELTA
         DELDIP(2,I)=(DELDIP(2,I)-DEL2(2))*0.5D0/DELTA
         DELDIP(3,I)=(DELDIP(3,I)-DEL2(3))*0.5D0/DELTA
         LL=LU+1
         LU=LL+I-1
         L=0
         IF(PRECIS)THEN
            DO 80 KOUNTF=LL,LU
               L=L+1
C
C       G2OLD = X + 1.0*DELTA
C       GROLD = X + 0.5*DELTA
C       GRAD  = X - 0.5*DELTA
C       G2RAD = X - 1.0*DELTA
C
               DUMY(L)= (8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
     1          /DELTA*FACT/24.D0
               EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
     1          /DELTA**3*FACT/56.D0
C
C  CORRECT FOR 4'TH ORDER CONTAMINATION
C
C#             CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
C#             DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
               FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
   80       CONTINUE
            L=L-1
            DO 90 K=I,NVAR
               L=L+1
               KK=(K*(K-1))/2+I
               DUMY(L)=(8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
     1          /DELTA*FACT/24.D0
               EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
     1          /DELTA**3*FACT/56.D0
C
C  CORRECT FOR 4'TH ORDER CONTAMINATION
C
C#             CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
C#             DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
               FMATRX(KK)=FMATRX(KK)+DUMY(L)
   90       CONTINUE
         ELSE
            DO 100 KOUNTF=LL,LU
               L=L+1
               DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
               FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
  100       CONTINUE
            L=L-1
            DO 110 K=I,NVAR
               L=L+1
               KK=(K*(K-1))/2+I
               DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
               FMATRX(KK)=FMATRX(KK)+DUMY(L)
  110       CONTINUE
         ENDIF
         IF(BIG)THEN
            WRITE(6,'(A)')' CONTRIBUTIONS TO F-MATRIX'
            WRITE(6,'(A)')' ELEMENT  +1.0*DELTA  +0.5*DELTA  -0.5*DEL'
     1//'TA  -1.0*DELTA   2''ND ORDER 4TH ORDER'
            WRITE(6,'(I7,6F12.6)')(L,G2OLD(L),GROLD(L),GRAD(L),G2RAD(L),
     1DUMY(L),EIGS(L),L=1,NVAR)
         ENDIF
         TIME3 = SECOND()
         TSTEP=TIME3-TIME2
         TLEFT= MAX(0.1D0,TLEFT-TSTEP)
         IF(TSTEP.GT.1.D6)TSTEP=TSTEP-1.D6
         TOTIME= TOTIME+TSTEP
         IF(RESFIL)THEN
            WRITE(6,'('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL ='
     1',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
         IF(LOG)WRITE(11,'('' STEP:'',I4,'' RESTART FILE WRITTEN, '',
     +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
            RESFIL=.FALSE.
         ELSE
            WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL =
     1'',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
        IF(LOG) WRITE(11,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, '',
     +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
         ENDIF
         ESTIM = TOTIME/I
         IF(TLAST-TLEFT.GT.TDUMP)THEN
            TLAST=TLEFT
            RESFIL=.TRUE.
            JSTART=1
            II=I
            CALL FORSAV(TOTIME,DELDIP,II,NVAR,FMATRX, COORD,NVAR,HEAT,
     1                EVECS,JSTART,FCONST)
         ENDIF
         IF(I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN
            WRITE(6,'(//10X,''- - - - - - - TIME UP - - - - - - -'',//)'
     1)
            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)
            WRITE(6,'(//10X,''FORCE MATRIX WRITTEN TO DISK'')')
            NREAL=-1
            RETURN
         ENDIF
  120 CONTINUE
      DO 130 I=1,NATOMS
         IF(ATMASS(I).LT.1.D-20.AND.LABELS(I).LT.99)THEN
            CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT,
     1                EVECS,ILOOP,FCONST)
            WRITE(6,'(A)')' AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART
     1'
            WRITE(6,'(A)')' FILE HAS BEEN WRITTEN AND THE JOB STOPPED'
            STOP
         ENDIF
  130 CONTINUE
      IF(ISTART.LE.NVAR .AND. INDEX(KEYWRD,'ISOT') .NE. 0)
     1CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT,
     2                EVECS,ILOOP,FCONST)
      RETURN
      END
