      SUBROUTINE DERIV(GEO,GRAD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'mopac.inc'
      DIMENSION GRAD(*), GEO(3,*)
      COMMON / EULER/ TVEC(3,3), ID
      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, DUMMY(MAXPAR)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /GRAVEC/ COSINE
      COMMON /GEOSYM/ NDEP, IDUMYS(MAXPAR,3)
      COMMON /PATH  / LATOM,LPARAM,REACT(200)
      COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
      COMMON /XYZGRA/ DXYZ(3,NUMATM*27)
      COMMON /ENUCLR/ ENUCLR
      COMMON /NUMCAL/ NUMCAL
      COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
      COMMON /WMATRX/ WJ(N2ELEC), WK(N2ELEC)
      COMMON /HMATRX/ H(MPACK)
      COMMON /ATHEAT/ ATHEAT
C***********************************************************************
C
C    DERIV CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
C          INTERNAL COORDINATES. THIS IS DONE BY FINITE DIFFERENCES.
C
C    THE MAIN ARRAYS IN DERIV ARE:
C        LOC    INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM
C               INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE
C               DERIVATIVE CALCULATION.
C        GEO    ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES.
C        GRAD   ON EXIT, CONTAINS THE DERIVATIVES
C
C***********************************************************************
      COMMON /KEYWRD / KEYWRD
      COMMON /ERRFN  / ERRFN(MAXPAR)
      CHARACTER*80 KEYWRD
      DIMENSION CHANGE(3), COORD(3,NUMATM), COLD(3,NUMATM*27)
     1,         XDERIV(3), XPARAM(MAXPAR), XJUC(3), W(N2ELEC)
      REAL WJ, WK
      LOGICAL DEBUG, TIMES, HALFE, FAST, SCF1, CI, PRECIS
      EQUIVALENCE (W,WJ)
      DATA ICALCN /0/
      IF(ICALCN.NE.NUMCAL) THEN
         I = INDEX(KEYWRD,'PRESS')
         PRESS=0.D0
         IF(I.NE.0) PRESS=READA(KEYWRD,I)*1476.8992D0
         IDLO=NATOMS+1
         IF(LABELS(NATOMS) .EQ. 107) THEN
            IDLO=NATOMS
            IF(LABELS(NATOMS-1) .EQ. 107)THEN
               IDLO=NATOMS-1
               IF(LABELS(NATOMS-2) .EQ. 107)THEN
                  IDLO=NATOMS-2
               ENDIF
            ENDIF
         ENDIF
         GRLIM=0.01D0
         DEBUG = (INDEX(KEYWRD,'DERIV') .NE. 0)
         PRECIS= (INDEX(KEYWRD,'PRECIS') .NE. 0)
         TIMES = (INDEX(KEYWRD,'TIME') .NE. 0)
         CI    = (INDEX(KEYWRD,'C.I.') .NE. 0)
         SCF1  = (INDEX(KEYWRD,'1SCF') .NE. 0)
         ICALCN=NUMCAL
         IF(INDEX(KEYWRD,'RESTART') .EQ. 0) THEN
            DO 10 I=1,NVAR
   10       ERRFN(I)=0.D0
         ENDIF
         GRLIM=0.01D0
         IF(PRECIS)GRLIM=0.0001D0
         IF(INDEX(KEYWRD,'FULSCF') .GT.0) GRLIM=1.D9
         HALFE = (NOPEN.GT.NCLOSE .OR. CI)
         IDELTA=-7
*
*   IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER
*
         IF(HALFE.AND.PRECIS) IDELTA=-3
         IF(HALFE.AND..NOT.PRECIS) IDELTA=-3
         FAST=.TRUE.
         CHANGE(1)= 10.D0**IDELTA
         CHANGE(2)= 10.D0**IDELTA
         CHANGE(3)= 10.D0**IDELTA
C
C    CHANGE(I) IS THE STEP SIZE USED IN CALCULATING THE DERIVATIVES.
C    FOR "CARTESIAN" DERIVATIVES, CALCULATED USING DCART,AN
C    INFINITESIMAL STEP, HERE 0.000001, IS ACCEPTABLE. IN THE
C    HALF-ELECTRON METHOD A QUITE LARGE STEP IS NEEDED AS FULL SCF
C    CALCULATIONS ARE NEEDED, AND THE DIFFERENCE BETWEEN THE TOTAL
C    ENERGIES IS USED. THE STEP CANNOT BE VERY LARGE, AS THE SECOND
C    DERIVITIVE IN FLEPO IS CALCULATED FROM THE DIFFERENCES OF TWO
C    FIRST DERIVATIVES. CHANGE(1) IS FOR CHANGE IN BOND LENGTH,
C    (2) FOR ANGLE, AND (3) FOR DIHEDRAL.
C
         XDERIV(1)= 0.5D0/CHANGE(1)
         XDERIV(2)= 0.5D0/CHANGE(2)
         XDERIV(3)= 0.5D0/CHANGE(3)
      ENDIF
      GNORM=0.D0
      IF(NVAR.EQ.0) RETURN
      IF(DEBUG)THEN
         WRITE(6,'('' GEO AT START OF DERIV'')')
         WRITE(6,'(F19.5,2F12.5)')((GEO(J,I),J=1,3),I=1,NATOMS)
      ENDIF
      DO 20 I=1,NVAR
         XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
   20 GNORM=GNORM+GRAD(I)**2
      GNORM=SQRT(GNORM)
      FAST=(GNORM .GT. GRLIM .AND. .NOT. SCF1 .OR. .NOT. HALFE)
      TIME1=ZECOND()
      IF(NDEP.NE.0) CALL SYMTRY
      CALL GMETRY(GEO,COORD)
      IF( .NOT. FAST ) THEN
         IF(DEBUG)WRITE(6,'('' DOING FULL SCF''''S IN DERIV'')')
         CALL HCORE(COORD,H,W, WJ, WK, ENUCLR)
         IF(NORBS*NELECS.GT.0)THEN
            CALL ITER(H, W, WJ, WK, AA,.TRUE.,.FALSE.)
         ELSE
            AA=0.D0
         ENDIF
         LINEAR=(NORBS*(NORBS+1))/2
         DO 30 I=1,LINEAR
   30    P(I)=PA(I)*2.D0
         AA=(AA+ENUCLR)
      ENDIF
      CALL DCART(COORD,DXYZ)
      IF(NDEP.NE.0) CALL SYMTRY
      CALL GMETRY(GEO,COORD)
      IJ=0
      DO 70 II=1,NUMAT
         DO 60 IL=L1L,L1U
            DO 60 JL=L2L,L2U
               DO 60 KL=L3L,L3U
                  DO 40 LL=1,3
   40             XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*JL+TVEC
     1(LL,3)*KL
                  IJ=IJ+1
                  DO 50 KK=1,3
                     COLD(KK,IJ)=XJUC(KK)
   50             CONTINUE
   60    CONTINUE
   70 CONTINUE
      SUM11=1.D-9
      SUM22=1.D-9
      SUM12=1.D-9
      DO 150 I=1,NVAR
         K=LOC(1,I)
         L=LOC(2,I)
         XSTORE=XPARAM(I)
         DO 80 J=1,NVAR
   80    GEO(LOC(2,J),LOC(1,J))=XPARAM(J)
         GEO(L,K)=XSTORE-CHANGE(L)
         IF(NDEP.NE.0) CALL SYMTRY
         CALL GMETRY(GEO,COORD)
C#         CALL GEOUT
C
C    USE LOOKUP TABLE OF CARTESIAN DERIVATIVES TO WORK OUT INTERNAL
C    COORDINATE DERIVATIVE.
C
         TOTL=0.D0
         IJ=0
         DO 130 II=1,NUMAT
            IF(ID.EQ.0) THEN
               DO 90 LL=1,3
   90          TOTL=TOTL+DXYZ(LL,II)*(COORD(LL,II)-COLD(LL,II))
            ELSE
               DO 120 IL=L1L,L1U
                  DO 120 JL=L2L,L2U
                     DO 120 KL=L3L,L3U
                        DO 100 LL=1,3
  100                   XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*J
     1L+TVEC(LL,3)*KL
                        IJ=IJ+1
                        DO 110 KK=1,3
                           TOTL=TOTL+DXYZ(KK,IJ)*(XJUC(KK)-COLD(KK,IJ))
  110                   CONTINUE
  120          CONTINUE
            ENDIF
  130    CONTINUE
         TOTL=TOTL*XDERIV(L)
C
C   IF NEEDED, CALCULATE "EXACT" DERIVITIVES.
C
         IF( .NOT. FAST ) THEN
            CALL HCORE(COORD,H,W, WJ, WK,ENUCLR)
            IF(NORBS*NELECS.GT.0)THEN
               CALL ITER(H,W, WJ, WK,EE,.TRUE.,.FALSE.)
            ELSE
               EE=0.D0
            ENDIF
            DO 140 II=1,LINEAR
  140       P(II)=PA(II)*2.D0
            EE=(EE+ENUCLR)
            TOTL1=(AA-EE)*23.061D0*XDERIV(L)*2.D0
C#            WRITE(6,*)AA-EE
            ERRFN(I)=TOTL1-TOTL
         ENDIF
         GEO(L,K)=XSTORE
         SUM11=SUM11+GRAD(I)**2
         SUM22=SUM22+TOTL**2
         SUM12=SUM12+TOTL*GRAD(I)
         GRAD(I)=TOTL+ERRFN(I)
  150 CONTINUE
      IF(DEBUG) THEN
         WRITE(6,'('' GRADIENTS'')')
         WRITE(6,'(10F8.3)')(GRAD(I),I=1,NVAR)
         WRITE(6,'('' ERROR FUNCTION'')')
         WRITE(6,'(10F8.3)')(ERRFN(I),I=1,NVAR)
      ENDIF
      COSINE=SUM12/SQRT(SUM11*SUM22)
      IF(DEBUG)
     1WRITE(6,'('' COSINE OF SEARCH DIRECTION ='',F30.6)')COSINE
      IF( .NOT. FAST ) COSINE=1.D0
      IF(TIMES)
     1WRITE(6,'('' TIME FOR DERIVATIVES'',F12.6)')ZECOND()-TIME1
      RETURN
      END
