      SUBROUTINE DERIV(GEO,GRAD)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      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, GETNAM*80
      dimension CHANGE(3), XDERIV(3), XJUC(3), AIDREF(MAXPAR)
      logical DEBUG, HALFE, SCF1, CI, PRECIS, SLOW, AIC, NOANCI,
     1       AIFRST, ISOK, GEOOK, INT
      save SCF1, HALFE, IDELTA, SLOW, XDERIV
      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=MIOA,file=GETNAM('FOR005'),
     1           status='OLD',BLANK='ZERO')
            rewind MIOA
C
C  ISOK IS SET FALSE: ONLY ONE SYSTEM ALLOWED
C
            ISOK=.FALSE.
            do 10 I=1,3
   10       read(MIOA,'(A)') LINE
            do 30 J=1,1000
               read(MIOA,'(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)
                  end if
   20          continue
************************************************************************
   30       if (INDEX(LINE,'AIDER').NE.0) go to 60
   40       write(MIOB,'(//,A)')' KEYWORD "AIDER" SPECIFIED, BUT NOT'
            write(MIOB,'(A)')' PRESENT AFTER Z-MATRIX.  JOB STOPPED'
            call filclo  ! J/M
            stop ' in DERIV - AIDER used, not pres after Z-mat.'
   50       write(MIOB,'(//,A)')
     %                   '  FAULT IN READ OF AB INITIO DERIVATIVES'
            write(MIOB,'(A)')'  DERIVATIVES READ IN ARE AS FOLLOWS'
            write(MIOB,'(6F12.6)') (AIDREF(J),J=1,I)
            call filclo  ! J/M
            stop ' in DERIV - fault in ab initio derivatives'
   60       continue
            if (NATOMS.GT.2) then
               J=3*NATOMS-6
            else
               J=1
            end if
            read(MIOA,*,end=50,ERR=50) (AIDREF(I),I=1,J)
            write(MIOB,'(/,A,/)')
     1' AB-INITIO DERIVATIVES IN KCAL/MOL/(ANGSTROM OR RADIAN)'
            write(MIOB,'(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
               else if (LOC(1,I).EQ.3) then
                  J=LOC(2,I)+1
               else
                  J=1
               end if
   70       AIDREF(I)=AIDREF(J)
            write(MIOB,'(/,A,/)')
     1' AB-INITIO DERIVATIVES FOR VARIABLES'
            write(MIOB,'(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.
     2                  LOC(2,I).EQ.3.AND.IDEPFN(J).EQ.14))
     3                AIDREF(I)=AIDREF(I)+SUM
   80             continue
   90          continue
               write(MIOB,'(/,A,/)')
     1 ' AB-INITIO DERIVATIVES AFTER SYMMETRY WEIGHTING'
               write(MIOB,'(5F12.6)') (AIDREF(J),J=1,NVAR)
            end if
         end if
         ICALCN=NUMCAL
         if (INDEX(KEYWRD,'RESTART') .EQ. 0) then
            do 100 I=1,NVAR
  100       ERRFN(I)=0.D0
         end if
         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)
      end if
      if (NVAR.EQ.0) return
      if (DEBUG) then
         write(MIOB,'('' GEO AT START OF DERIV'')')
         write(MIOB,'(F19.5,2F12.5)') ((GEO(J,I),J=1,3),I=1,NATOMS)
      end if
      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))
      end if
      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(MIOB,*) 'DOING ANALYTICAL C.I. DERIVATIVES'
         call DERNVO(COORD,DXYZ)
      else
         if (DEBUG)
     %     write(MIOB,*) 'DOING VARIATIONALLY OPIMIZED DERIVATIVES'
         call DCART(COORD,DXYZ)
      end if
      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
      end if
      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.
     1       DOT(DXYZ,DXYZ,3*NUMAT).GT.MAX(4.D0,SUM*4.D0)) 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(MIOB,'(//,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(MIOB,'(//,A)')'     CURRENT GEOMETRY'
                  call GEOUT(MIOB)
                  call filclo  ! J/M
                  stop ' in DERIV - cartesian vs. local coords problem'
               end if
  170       continue
         end if
      end if
C
C  THIS CODE IS ONLY USED if THE KEYWORD NOANCI IS SPECIFIED
      if (SLOW) then
         if (DEBUG) write(MIOB,*) '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)
      end if
      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)
         end if
C#         write(MIOB,'('' GRADIENTS BEFORE AI CORRECTION'')')
C#         write(MIOB,'(10F8.3)') (GRAD(I),I=1,NVAR)
         do 210 I=1,NVAR
  210    GRAD(I)=GRAD(I)+AICORR(I)
      end if
  220 if (DEBUG) then
         write(MIOB,'('' GRADIENTS'')')
         write(MIOB,'(10F8.3)') (GRAD(I),I=1,NVAR)
         if (SLOW) then
            write(MIOB,'('' ERROR FUNCTION'')')
            write(MIOB,'(10F8.3)') (ERRFN(I),I=1,NVAR)
         end if
      end if
      if (DEBUG)
     1  write(MIOB,'('' COSINE OF SEARCH DIRECTION ='',F30.6)')
     2        COSINE
      return
      end
