      SUBROUTINE DERIV(GEO,GRAD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION GRAD(*), GEO(3,*)
************************************************************************
*
*    DERIV CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
*          INTERNAL COORDINATES. THIS IS DONE BY FINITE DIFFERENCES.
*
*    THE MAIN ARRAYS IN DERIV ARE:
*        LOC    INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM
*               INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE
*               DERIVATIVE CALCULATION.
*        GEO    ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES.
*        GRAD   ON EXIT, CONTAINS THE DERIVATIVES
*
************************************************************************
      COMMON / EULER/ TVEC(3,3), ID
      COMMON /OKMANY/ ISOK
      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, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1LOCDEP(MAXPAR)
      COMMON /PATH  / LATOM,LPARAM,REACT(200)
      COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
      COMMON /XYZGRA/ DXYZ(9*NUMATM)
      COMMON /ENUCLR/ ENUCLR
      COMMON /NUMCAL/ NUMCAL
      COMMON /HMATRX/ H(MPACK)
      COMMON /ATHEAT/ ATHEAT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ERRFN / ERRFN(MAXPAR), AICORR(MAXPAR)
      COMMON /WORK1 / WORK2(22*MPACK)
      COMMON /GENRAL/ COORD(3,NUMATM), COLD(3,NUMATM*3), GOLD(MAXPAR),
     1 XPARAM(MAXPAR)
      CHARACTER*241 KEYWRD, LINE*80
      DIMENSION CHANGE(3), XDERIV(3), XJUC(3), AIDREF(MAXPAR)
      SAVE SCF1, HALFE, IDELTA, SLOW, XDERIV
      LOGICAL DEBUG, HALFE, SCF1, CI, PRECIS, SLOW, AIC, NOANCI,
     1AIFRST, ISOK, GEOOK, INT
      DATA ICALCN /0/
      IF(ICALCN.NE.NUMCAL) THEN
         AIFRST= (INDEX(KEYWRD,'RESTART').EQ.0)
         DEBUG = (INDEX(KEYWRD,'DERIV') .NE. 0)
         PRECIS= (INDEX(KEYWRD,'PREC') .NE. 0)
         INT   = (INDEX(KEYWRD,' XYZ') .EQ. 0)
         GEOOK = (INDEX(KEYWRD,'GEO-OK') .NE. 0)
         CI    = (INDEX(KEYWRD,'C.I.') .NE. 0)
         SCF1  = (INDEX(KEYWRD,'1SCF') .NE. 0)
         AIC=(INDEX(KEYWRD,'AIDER').NE.0)
         ICAPA=ICHAR('A')
         ILOWA=ICHAR('a')
         ILOWZ=ICHAR('z')
         IF(AIC.AND.AIFRST)THEN
            OPEN(UNIT=5,FILE='FOR005',STATUS='OLD',BLANK='ZERO')
            REWIND 5
C
C  ISOK IS SET FALSE: ONLY ONE SYSTEM ALLOWED
C
            ISOK=.FALSE.
            DO 10 I=1,3
   10       READ(5,'(A)')LINE
            DO 30 J=1,1000
               READ(5,'(A)',END=40,ERR=40)LINE
************************************************************************
               DO 20 I=1,80
                  ILINE=ICHAR(LINE(I:I))
                  IF(ILINE.GE.ILOWA.AND.ILINE.LE.ILOWZ) THEN
                     LINE(I:I)=CHAR(ILINE+ICAPA-ILOWA)
                  ENDIF
   20          CONTINUE
************************************************************************
   30       IF(INDEX(LINE,'AIDER').NE.0)GOTO 60
   40       WRITE(6,'(//,A)')' KEYWORD "AIDER" SPECIFIED, BUT NOT'
            WRITE(6,'(A)')' PRESENT AFTER Z-MATRIX.  JOB STOPPED'
            STOP
   50       WRITE(6,'(//,A)')'  FAULT IN READ OF AB INITIO DERIVATIVES'
            WRITE(6,'(A)')'  DERIVATIVES READ IN ARE AS FOLLOWS'
            WRITE(6,'(6F12.6)')(AIDREF(J),J=1,I)
            STOP
   60       CONTINUE
            IF(NATOMS.GT.2)THEN
               J=3*NATOMS-6
            ELSE
               J=1
            ENDIF
            READ(5,*,END=50,ERR=50)(AIDREF(I),I=1,J)
            WRITE(6,'(/,A,/)')
     1' AB-INITIO DERIVATIVES IN KCAL/MOL/(ANGSTROM OR RADIAN)'
            WRITE(6,'(5F12.6)')(AIDREF(I),I=1,J)
            DO 70 I=1,NVAR
               IF(LOC(1,I).GT.3)THEN
                  J=3*LOC(1,I)+LOC(2,I)-9
               ELSEIF(LOC(1,I).EQ.3)THEN
                  J=LOC(2,I)+1
               ELSE
                  J=1
               ENDIF
   70       AIDREF(I)=AIDREF(J)
            WRITE(6,'(/,A,/)')
     1' AB-INITIO DERIVATIVES FOR VARIABLES'
            WRITE(6,'(5F12.6)')(AIDREF(I),I=1,NVAR)
            IF(NDEP.NE.0)THEN
               DO 90 I=1,NVAR
                  SUM=AIDREF(I)
                  DO 80 J=1,NDEP
                     IF(LOC(1,I).EQ.LOCPAR(J).AND.(LOC(2,I).EQ.IDEPFN(J)
     1.OR.LOC(2,I).EQ.3.AND.IDEPFN(J).EQ.14)) AIDREF(I)=AIDREF(I)+SUM
   80             CONTINUE
   90          CONTINUE
               WRITE(6,'(/,A,/)')
     1' AB-INITIO DERIVATIVES AFTER SYMMETRY WEIGHTING'
               WRITE(6,'(5F12.6)')(AIDREF(J),J=1,NVAR)
            ENDIF
         ENDIF
         ICALCN=NUMCAL
         IF(INDEX(KEYWRD,'RESTART') .EQ. 0) THEN
            DO 100 I=1,NVAR
  100       ERRFN(I)=0.D0
         ENDIF
         GRLIM=0.01D0
         IF(PRECIS)GRLIM=0.0001D0
         HALFE = (NOPEN.GT.NCLOSE.AND.FRACT.NE.2.D0.AND.FRACT.NE.0.D0
     1 .OR. CI)
         IDELTA=-7
*
*   IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER
*
         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
      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
      GNORM=0.D0
      DO 110 I=1,NVAR
         GOLD(I)=GRAD(I)
         XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
  110 GNORM=GNORM+GRAD(I)**2
      GNORM=SQRT(GNORM)
      SLOW=.FALSE.
      NOANCI=.FALSE.
      IF(HALFE) THEN
         NOANCI=(INDEX(KEYWRD,'NOANCI').NE.0 .OR. NOPEN.EQ.NORBS)
         SLOW=(NOANCI.AND.(GNORM .LT. GRLIM .OR. SCF1))
      ENDIF
      IF(NDEP.NE.0) CALL SYMTRY
      CALL GMETRY(GEO,COORD)
C
C  COORD NOW HOLDS THE CARTESIAN COORDINATES
C
      IF(HALFE.AND..NOT.NOANCI) THEN
         IF(DEBUG)WRITE(6,*) 'DOING ANALYTICAL C.I. DERIVATIVES'
         CALL DERNVO(COORD,DXYZ)
      ELSE
         IF(DEBUG)WRITE(6,*) 'DOING VARIATIONALLY OPIMIZED DERIVATIVES'
         CALL DCART(COORD,DXYZ)
      ENDIF
      IJ=0
      DO 150 II=1,NUMAT
         DO 140 IL=L1L,L1U
            DO 140 JL=L2L,L2U
               DO 140 KL=L3L,L3U
C$DOIT ASIS
                  DO 120 LL=1,3
  120             XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*JL+TVEC
     1(LL,3)*KL
                  IJ=IJ+1
C$DOIT ASIS
                  DO 130 KK=1,3
                     COLD(KK,IJ)=XJUC(KK)
  130             CONTINUE
  140    CONTINUE
  150 CONTINUE
      STEP=CHANGE(1)
      CALL JCARIN (COORD,XPARAM,STEP,PRECIS,WORK2,NCOL)
      CALL MXM (WORK2,NVAR,DXYZ,NCOL,GRAD,1)
      IF (PRECIS) THEN
         STEP=0.5D0/STEP
      ELSE
         STEP=1.0D0/STEP
      ENDIF
      DO 160 I=1,NVAR
  160 GRAD(I)=GRAD(I)*STEP
C
C  NOW TO ENSURE THAT INTERNAL DERIVATIVES ACCURATELY REFLECT CARTESIAN
C  DERIVATIVES
C
      IF(INT.AND. .NOT. GEOOK .AND. NVAR.GE.NUMAT*3-6.AND.ID.EQ.0)THEN
C
C  NUMBER OF VARIABLES LOOKS O.K.
C
         SUM=DOT(GRAD,GRAD,NVAR)
         IF(SUM.LT.2.D0.AND.DOT(DXYZ,DXYZ,3*NUMAT).GT.MAX(4.D0,SUM*4.D0)
     1)THEN
C
C OOPS, LOOKS LIKE AN ERROR.
C
            DO 170 I=1,NVAR
               J=XPARAM(I)/3.141D0
               IF(LOC(2,I).EQ.2.AND.LOC(1,I).GT.3.AND.
     1 ABS(XPARAM(I)-J*3.1415926D0).LT.0.005D0 )THEN
C
C  ERROR LOCATED, BUT CANNOT CORRECT IN THIS RUN
C
                  WRITE(6,'(//,3(A,/),I3,A)')
     1' INTERNAL COORDINATE DERIVATIVES DO NOT REFLECT',
     2' CARTESIAN COORDINATE DERIVATIVES',
     3' TO CORRECT ERROR, INCREASE DIHEDRAL OF ATOM',LOC(1,I),
     4' BY 90 DEGREES'
                  WRITE(6,'(//,A)')'     CURRENT GEOMETRY'
                  CALL GEOUT(6)
                  STOP
               ENDIF
  170       CONTINUE
         ENDIF
      ENDIF
C
C  THIS CODE IS ONLY USED IF THE KEYWORD NOANCI IS SPECIFIED
      IF(SLOW)THEN
         IF(DEBUG)WRITE(6,*) 'DOING FULL SCF DERIVATIVES'
         CALL DERITR(ERRFN,GEO)
C
C THE ARRAY ERRFN HOLDS THE EXACT DERIVATIVES MINUS THE APPROXIMATE
C DERIVATIVES
         DO 180 I=1,NVAR
  180    ERRFN(I)=ERRFN(I)-GRAD(I)
      ENDIF
      COSINE=DOT(GRAD,GOLD,NVAR)/
     1SQRT(DOT(GRAD,GRAD,NVAR)*DOT(GOLD,GOLD,NVAR)+1.D-20)
      DO 190 I=1,NVAR
  190 GRAD(I)=GRAD(I)+ERRFN(I)
      IF(AIC)THEN
         IF(AIFRST)THEN
            AIFRST=.FALSE.
            DO 200 I=1,NVAR
  200       AICORR(I)=-AIDREF(I)-GRAD(I)
         ENDIF
C#         WRITE(6,'('' GRADIENTS BEFORE AI CORRECTION'')')
C#         WRITE(6,'(10F8.3)')(GRAD(I),I=1,NVAR)
         DO 210 I=1,NVAR
  210    GRAD(I)=GRAD(I)+AICORR(I)
      ENDIF
  220 IF(DEBUG) THEN
         WRITE(6,'('' GRADIENTS'')')
         WRITE(6,'(10F8.3)')(GRAD(I),I=1,NVAR)
         IF(SLOW)THEN
            WRITE(6,'('' ERROR FUNCTION'')')
            WRITE(6,'(10F8.3)')(ERRFN(I),I=1,NVAR)
         ENDIF
      ENDIF
      IF(DEBUG)
     1WRITE(6,'('' COSINE OF SEARCH DIRECTION ='',F30.6)')COSINE
      RETURN
      END
