      SUBROUTINE FMAT(FMATRX, NREAL, TSCF, TDER, DELDIP, HEAT)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      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),
     1 GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3), G2OLD(MAXPAR)
     2, EIGS(MAXPAR), G2RAD(MAXPAR),
     3 FCONST(MAXPAR)
      character*241 KEYWRD
      logical DEBUG, RESTRT, PRNT, RESFIL, ANALYT, PRECIS, BIG, LOG
      character SPACE*1, DOTT*1, ZERO*1, NINE*1
      save  SPACE, DOTT, ZERO, NINE, FACT
      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)
         end if
   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(MIOB,'(//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
      end if
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
      end if
      if (TSCF.GT.0)
     1   write(MIOB,'(/10X,''ESTIMATED TIME TO COMPLETE CALCULATION ='',
     2                F9.2,'' SECONDS'')')  ESTIME
      if (RESTRT) then
         if (ISTART.LE.NVAR)
     1    write(MIOB,'(/10X,''STARTING AGAIN AT LINE'',18X,I4)')
     2                     ISTART
         write(MIOB,'(/10X,''TIME USED UP TO RESTART ='',F22.2)')
     %                     TOTIME
      end if
      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(MIOB,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELTA
            G2OLD(1)=100.D0
            COORDL(I)=COORDL(I)+DELTA
            call COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
            if (DEBUG)  write(MIOB,'(A,F12.5)')  ' GNORM +1.0*DELTA',
     %                                   SQRT(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(MIOB,'(A,F12.5)') ' GNORM -1.0*DELTA',
     %                                   SQRT(DOT(G2RAD,G2RAD,NVAR))
         else
            if (DEBUG)
     %        write(MIOB,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELTA
         end if
         COORDL(I)=COORDL(I)+0.5D0*DELTA
         GROLD(1)=100.D0
         call COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.)
         if (DEBUG)  write(MIOB,'(A,F12.5)')  ' GNORM +0.5*DELTA',
     %                                   SQRT(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(MIOB,'(A,F12.5)')  ' GNORM -0.5*DELTA',
     %                                   SQRT(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
         end if
         if (BIG) then
            write(MIOB,'(A)') ' CONTRIBUTIONS TO F-MATRIX'
            write(MIOB,'(A)')
     %  ' ELEMENT  +1.0*DELTA  +0.5*DELTA  -0.5*DEL'//
     %  'TA  -1.0*DELTA   2''ND ORDER 4TH ORDER'
            write(MIOB,'(I7,6F12.6)')
     %                  (L,G2OLD(L),GROLD(L),GRAD(L),G2RAD(L),
     %                                        DUMY(L),EIGS(L),L=1,NVAR)
         end if
         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(MIOB,
     %       '('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL ='',
     %   F10.2,'' TIME LEFT:'',F10.2)')
     %                                  I,TOTIME,TLEFT
         if (LOG) write(MIOE,
     %    '('' STEP:'',I4,'' RESTART FILE WRITTEN, '',
     %     ''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')
     %     I,TOTIME,TLEFT
            RESFIL=.FALSE.
         else
            write(MIOB,
     %       '('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL ='',
     %   F10.2,'' TIME LEFT:'',F10.2)')
     %                                  I,TSTEP,TOTIME,TLEFT
        if (LOG) write(MIOE,'('' STEP:'',I4,'' TIME ='',F9.2,
     %    '' SECS, '',''INTEGRAL ='',F10.2,
     %    '' TIME LEFT:'',F10.2)')
     %                              I,TSTEP,TOTIME,TLEFT
         end if
         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)
         end if
         if (I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) then
            write(MIOB,
     %'(//10X,''- - - - - - - TIME UP - - - - - - -'',//)')
            write(MIOB,'(/10X,'' POINT REACHED ='',I4)') I
            write(MIOB,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')')
            write(MIOB,
     %         '(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2,
     %         '' SECONDS'')')    ESTIM
            JSTART=1
            II=I
            call FORSAV(TOTIME,DELDIP,II,NVAR,FMATRX, COORD,NVAR,HEAT,
     1                EVECS,JSTART,FCONST)
            write(MIOB,'(//10X,''FORCE MATRIX WRITTEN TO DISK'')')
            NREAL=-1
            return
         end if
  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(MIOB,'(A)')
     %                 ' AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART'
            write(MIOB,'(A)')
     %                 ' FILE HAS BEEN WRITTEN AND THE JOB STOPPED'
            call filclo  !  J/M
            stop ' in FMAT - A zero mass atom found!'
         end if
  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
