SUBROUTINE DERI1(C,NORBS,COORD,NUMBER,WORK,GRAD
     1                 ,F,MINEAR,FD,WMAT,HMAT,FMAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION WMAT(MPACK),HMAT(MPACK*2),FMAT(MPACK*2)
*********************************************************************
*
*     DERI1 COMPUTE THE NON-RELAXED DERIVATIVE OF THE NON-VARIATIONALLY
*     OPTIMIZED WAVEFUNCTION ENERGY WITH RESPECT TO ONE CARTESIAN
*     COORDINATE AT A TIME
*                             AND
*     COMPUTE THE NON-RELAXED FOCK MATRIX DERIVATIVE IN M.O BASIS AS
*     REQUIRED IN THE RELAXATION SECTION (ROUTINE 'DERI2').
*
*   INPUT
*     C(NORBS,NORBS) : M.O. COEFFICIENTS.
*     COORD  : CARTESIAN COORDINATES ARRAY.
*     NUMBER : LOCATION OF THE REQUIRED VARIABLE IN COORD.
*     WORK   : WORK ARRAY OF SIZE N*N.
*     WMAT     : WORK ARRAYS FOR d<PQ|RS> (2-CENTERS  A.O)
*   OUTPUT
*     C,COORD,NUMBER : NOT MODIFIED.
*     GRAD   : DERIVATIVE OF THE HEAT OF FORMATION WITH RESPECT TO
*              COORD(NUMBER), WITHOUT RELAXATION CORRECTION.
*     F(MINEAR) : NON-RELAXED FOCK MATRIX DERIVATIVE WITH RESPECT TO
*              COORD(NUMBER), EXPRESSED IN M.O BASIS, SCALED AND
*              PACKED, OFF-DIAGONAL BLOCKS ONLY.
*     FD     : IDEM BUT UNSCALED, DIAGONAL BLOCKS, C.I-ACTIVE ONLY.
*
************************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM)
     1               ,NLAST(NUMATM), NDUMY1, NELECS,NALPHA,NBETA
     2               ,NCLOSE,NOPEN,NDUMY,FRACT
     3       /VECTOR/ CDUM(MORB2),EIGS(MAXORB),WDUM(MORB2),EIGB(MAXORB)
     4       /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
      COMMON /CIBITS/ NMOS,LAB,NELEC,NBO(3)
     1       /HMATRX/ H(MPACK)
     2       /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
     3       /CIVECT/ CONF(NMECI**4+NMECI**2)
      COMMON /FOKMAT/ FDUMY(MPACK), SCALAR(MPACK)
      COMMON /NVOMAT/ DIAG(MPACK/2)
     1       /KEYWRD/ KEYWRD
      COMMON /NUMCAL/ NUMCAL
      DIMENSION COORD(*),C(NORBS,NORBS),WORK(NORBS,NORBS),F(*),FD(*)
      CHARACTER KEYWRD*241
      LOGICAL DEBUG
      DATA ICALCN /0/
C
      IF(ICALCN.NE.NUMCAL) THEN
         DEBUG=INDEX(KEYWRD,'DERI1').NE.0
         IPRT=6
         LINEAR=NORBS*(NORBS+1)/2
         ICALCN=NUMCAL
      ENDIF
      IF(DEBUG) CALL TIMER('BEFORE DERI1')
      STEP=1.D-3
C
C     2 POINTS FINITE DIFFERENCE TO GET THE INTEGRAL DERIVATIVES
C     ----------------------------------------------------------
C     STORED IN HMAT AND WMAT, WITHOUT DIVIDING BY THE STEP SIZE.
C
      NATI=(NUMBER-1)/3+1
      NATX=NUMBER-3*(NATI-1)
      CALL DHCORE (COORD,HMAT,WMAT,ENUCL2,NATI,NATX,STEP)
C
C HMAT HOLDS THE ONE-ELECTRON DERIVATIVES OF ATOM NATI FOR DIRECTION
C      NATX W.R.T. ALL OTHER ATOMS
C WMAT HOLDS THE TWO-ELECTRON DERIVATIVES OF ATOM NATI FOR DIRECTION
C      NATX W.R.T. ALL OTHER ATOMS
      STEP=0.5D0/STEP
C
C     NON-RELAXED FOCK MATRIX DERIVATIVE IN A.O BASIS.
C     ------------------------------------------------
C     STORED IN FMAT, DIVIDED BY STEP.
C
      CALL SCOPY (LINEAR,HMAT,1,FMAT,1)
      CALL DFOCK2 (FMAT,P,PA,WMAT,NUMAT,NFIRST,NMIDLE,NLAST,NATI)
C
C  FMAT HOLDS THE ONE PLUS TWO - ELECTRON DERIVATIVES OF ATOM NATI FOR
C       DIRECTION NATX W.R.T. ALL OTHER ATOMS
C
C       DERIVATIVE OF THE SCF-ONLY ENERGY (I.E BEFORE C.I CORRECTION)
C
      GRAD=(HELECT(NORBS,P,HMAT,FMAT)+ENUCL2)*STEP
C     TAKE STEP INTO ACCOUNT IN FMAT
      DO 10 I=1,LINEAR
   10 FMAT(I)=FMAT(I)*STEP
C
C     RIGHT-HAND SIDE SUPER-VECTOR F = C' FMAT C USED IN RELAXATION
C     -----------------------------------------------------------
C     STORED IN NON-STANDARD PACKED FORM IN F(MINEAR) AND FD.
C     THE SUPERVECTOR IS THE NON-RELAXED FOCK MATRIX DERIVATIVE IN
C     M.O BASIS: F(IJ)= ( (C' * FOCK * C)(I,J) )   WITH I.GT.J .
C     F IS SCALED AND PACKED IN SUPERVECTOR FORM WITH
C                THE CONSECUTIVE FOLLOWING OFF-DIAGONAL BLOCKS:
C             1) OPEN-CLOSED  I.E. F(IJ)=F(I,J) WITH I OPEN & J CLOSED
C                                  AND I RUNNING FASTER THAN J,
C             2) VIRTUAL-CLOSED SAME RULE OF ORDERING,
C             3) VIRTUAL-OPEN   SAME RULE OF ORDERING.
C     FD IS PACKED OVER THE C.I-ACTIVE M.O WITH
C                THE CONSECUTIVE DIAGONAL BLOCKS:
C             1) CLOSED-CLOSED   IN CANONICAL ORDER, WITHOUT THE
C                                DIAGONAL ELEMENTS,
C             2) OPEN-OPEN       SAME RULE OF ORDERING,
C             3) VIRTUAL-VIRTUAL SAME RULE OF ORDERING.
C
C     PART 1 : WORK(N,N) = FMAT(N,N) * C(N,N)
      DO 20 I=1,NORBS
   20 CALL SUPDOT (WORK(1,I),FMAT,C(1,I),NORBS,1)
C
C     PART 2 : F(IJ) =  (C' * WORK)(I,J) ... OFF-DIAGONAL BLOCKS.
      L=1
      IF(NBO(2).NE.0 .AND. NBO(1).NE.0) THEN
C        OPEN-CLOSED
         CALL MTXM (C(1,NBO(1)+1),NBO(2),WORK,NORBS,F(L),NBO(1))
         L=L+NBO(2)*NBO(1)
      ENDIF
      IF(NBO(3).NE.0 .AND. NBO(1).NE.0) THEN
C        VIRTUAL-CLOSED
         CALL MTXM (C(1,NOPEN+1),NBO(3),WORK,NORBS,F(L),NBO(1))
         L=L+NBO(3)*NBO(1)
      ENDIF
      IF(NBO(3).NE.0 .AND. NBO(2).NE.0) THEN
C        VIRTUAL-OPEN
         CALL MTXM (C(1,NOPEN+1),NBO(3),WORK(1,NBO(1)+1),NORBS,F(L),NBO(
     12))
      ENDIF
C     SCALE F ACCORDING TO THE DIAGONAL METRIC TENSOR 'SCALAR '.
      DO 30 I=1,MINEAR
   30 F(I)=F(I)*SCALAR(I)
      IF(DEBUG)THEN
         WRITE(6,*)' F IN DERI1'
         J=MIN(20,MINEAR)
         WRITE(6,'(5F12.6)')(F(I),I=1,J)
      ENDIF
C
C     PART 3 : SUPER-VECTOR FD, C.I-ACTIVE DIAGONAL BLOCKS, UNSCALED.
      L=1
      NEND=0
      DO 50 LOOP=1,3
         NINIT=NEND+1
         NEND =NEND+NBO(LOOP)
         N1=MAX(NINIT,NELEC+1   )
         N2=MIN(NEND ,NELEC+NMOS)
         IF(N2.LT.N1) GO TO 50
         DO 40 I=N1,N2
            IF(I.GT.NINIT) THEN
               CALL MXM (C(1,I),1,WORK(1,NINIT),NORBS,FD(L),I-NINIT)
               L=L+I-NINIT
            ENDIF
   40    CONTINUE
   50 CONTINUE
C
C     NON-RELAXED C.I CORRECTION TO THE ENERGY DERIVATIVE.
C     ----------------------------------------------------
C
C     C.I-ACTIVE FOCK EIGENVALUES DERIVATIVES, STORED IN FD(CONTINUED).
      LCUT=L
      DO 60 I=NELEC+1,NELEC+NMOS
         FD(L)=DOT(C(1,I),WORK(1,I),NORBS)
   60 L=L+1
C
C     C.I-ACTIVE 2-ELECTRONS INTEGRALS DERIVATIVES. STORED IN XY.
C   FMAT IS USED HERE AS SCRATCH SPACE
C
      CALL DIJKL1 (C(1,NELEC+1),NORBS,NATI,WMAT,FMAT,HMAT,FMAT)
      DO 70 I=1,NMOS
         DO 70 J=1,NMOS
            DO 70 K=1,NMOS
               DO 70 L=1,NMOS
   70 XY(I,J,K,L)=XY(I,J,K,L)*STEP
C
C     BUILD THE C.I MATRIX DERIVATIVE, STORED IN WMAT.
      CALL MECID (FD(LCUT-NELEC),GSE,EIGB,WORK)
      IF(DEBUG)THEN
         WRITE(6,*)' GSE:',GSE
C#      WRITE(6,*)' EIGB:',(EIGB(I),I=1,10)
C#      WRITE(6,*)' WORK:',(WORK(I,1),I=1,10)
      ENDIF
      CALL MECIH (WORK,WMAT,NMOS,LAB)
C
C     NON-RELAXED C.I CONTRIBUTION TO THE ENERGY DERIVATIVE.
      CALL SUPDOT (WORK,WMAT,CONF,LAB,1)
      GRAD=(GRAD+DOT(CONF,WORK,LAB))*23.061D0
      IF(DEBUG) THEN
         WRITE(IPRT,'('' * * * GRADIENT COMPONENT NUMBER'',I4)')NUMBER
         WRITE(IPRT,'('' NON-RELAXED C.I-ACTIVE FOCK EIGENVALUES '',
     1                ''DERIVATIVES (E.V.)'')')
         WRITE(IPRT,'(8F10.4)')(FD(LCUT-1+I),I=1,NMOS)
         WRITE(IPRT,'('' NON-RELAXED 2-ELECTRONS DERIVATIVES (E.V.)''/
     1''    I    J    K    L       d<I(1)J(1)|K(2)L(2)>'')')
         DO 80 I=1,NMOS
            DO 80 J=1,I
               DO 80 K=1,I
                  LL=K
                  IF(K.EQ.I) LL=J
                  DO 80 L=1,LL
   80    WRITE(IPRT,'(4I5,F20.10)')
     1              NELEC+I,NELEC+J,NELEC+K,NELEC+L,XY(I,J,K,L)
         WRITE(IPRT,'('' NON-RELAXED GRADIENT COMPONENT'',F10.4,
     1'' KCAL/MOLE'')')GRAD
         CALL TIMER('AFTER DERI1')
      ENDIF
      RETURN
      END
      SUBROUTINE DHCORE (COORD,H,W,ENUCLR,NATI,NATX,STEP)
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION COORD(3,*),H(*),W(*)
C
C  DHCORE GENERATES THE 1-ELECTRON  AND 2-ELECTRON INTEGRALS DERIVATIVES
C         WITH RESPECT TO THE CARTESIAN COORDINATE COORD (NATX,NATI).
C
C  INPUT
C      COORD     : CARTESIAN  COORDINATES OF THE MOLECULE.
C      NATI,NATX : INDICES OF THE MOVING COORDINATE.
C      STEP      : STEP SIZE OF THE 2-POINTS FINITE DIFFERENCE.
C  OUTPUT
C      H         : 1-ELECTRON INTEGRALS DERIVATIVES (PACKED CANONICAL).
C      W         : 2-ELECTRON INTEGRALS DERIVATIVES (ORDERED AS REQUIRED
C                             IN DFOCK2 AND DIJKL1).
C      ENUCLR    : NUCLEAR ENERGY DERIVATIVE.
C
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
     3       /MOLORB/ USPD(MAXORB),DUMY(MAXORB)
     4       /KEYWRD/ KEYWRD
     5       /WMATRX/ WDUMMY(N2ELEC*2)
      CHARACTER*241 KEYWRD
      LOGICAL FIRST,MINDO
      DIMENSION E1B(10),DE1B(10),E2A(10),DE2A(10)
     1         ,DI(9,9),DDI(9,9),WJD(101),DWJD(101)
     2,NB(0:8)
      DATA NB/1,0,0,10,0,0,0,0,45/
      DATA FIRST/.TRUE./
      IF (FIRST) THEN
         IONE=1
         CUTOFF=1.D10
         FIRST=.FALSE.
         MINDO=INDEX(KEYWRD,'MINDO') .NE. 0
      ENDIF
      DO 10 I=1,(NORBS*(NORBS+1))/2
   10 H(I)=0
      ENUCLR=0.D0
      KR=1
      NROW=0
      I=NATI
      CSAVE=COORD(NATX,NATI)
      IA=NFIRST(NATI)
      IB=NLAST(NATI)
      IC=NMIDLE(NATI)
      NI=NAT(NATI)
      NROW=-NB(IB-IA)
      DO 20 J=1,NUMAT
   20 NROW=NROW+NB(NLAST(J)-NFIRST(J))
      NCOL=NB(NLAST(NATI)-NFIRST(NATI))
      NBAND2=0
      DO 120 J=1,NUMAT
         IF (J.EQ.NATI) GO TO 120
         JA=NFIRST(J)
         JB=NLAST(J)
         JC=NMIDLE(J)
         NJ=NAT(J)
         COORD(NATX,NATI)=CSAVE+STEP
         CALL H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DI)
C
C  THE FOLLOWING STYLE WAS NECESSARY TO GET ROUND A BUG IN THE
C  GOULD COMPILER
C
         COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
         CALL H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DDI)
C
C     FILL THE ATOM-OTHER ATOM ONE-ELECTRON MATRIX.
C
         I2=0
         IF (IA.GT.JA) THEN
            DO 30 I1=IA,IB
               IJ=I1*(I1-1)/2+JA-1
               I2=I2+1
               J2=0
               DO 30 J1=JA,JB
                  IJ=IJ+1
                  J2=J2+1
   30       H(IJ)=H(IJ)+(DI(I2,J2)-DDI(I2,J2))
         ELSE
            DO 40 I1=JA,JB
               IJ=I1*(I1-1)/2+IA-1
               I2=I2+1
               J2=0
               DO 40 J1=IA,IB
                  IJ=IJ+1
                  J2=J2+1
   40       H(IJ)=H(IJ)+(DI(J2,I2)-DDI(J2,I2))
         ENDIF
C
C     CALCULATE THE TWO-ELECTRON INTEGRALS, W; THE ELECTRON NUCLEAR TERM
C     E1B AND E2A; AND THE NUCLEAR-NUCLEAR TERM ENUC.
C
         KRO=KR
         NBAND1=NBAND2+1
         NBAND2=NBAND2+NB(NLAST(J)-NFIRST(J))
         IF (MINDO) THEN
            COORD(NATX,NATI)=CSAVE+STEP
            CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
            KR=KRO
            COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
            CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
            IF (KR.GT.KRO) THEN
               DO 50 K=1,KR-KRO+1
   50          W(KRO+K-1)=WJD(K)-DWJD(K)
            ENDIF
         ELSE
            COORD(NATX,NATI)=CSAVE+STEP
            CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
            KR=KRO
            COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
            CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
            IF (KR.GT.KRO) THEN
               DO 60 K=1,KR-KRO+1
   60          WJD(K)=WJD(K)-DWJD(K)
               J7=0
               DO 70 I1=KRO,KR
                  J7=J7+1
   70          W(I1)=WJD(J7)
            ENDIF
         ENDIF
         COORD(NATX,NATI)=CSAVE
         ENUCLR = ENUCLR + ENUC-DENUC
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM I.
C
         I2=0
         DO 80 I1=IA,IC
            II=I1*(I1-1)/2+IA-1
            DO 80 J1=IA,I1
               II=II+1
               I2=I2+1
   80    H(II)=H(II)+E1B(I2)-DE1B(I2)
C     CONTRIB D, CNDO.
         DO 90 I1=IC+1,IB
            II=(I1*(I1+1))/2
   90    H(II)=H(II)+E1B(1)-DE1B(1)
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM J.
C
         I2=0
         DO 100 I1=JA,JC
            II=I1*(I1-1)/2+JA-1
            DO 100 J1=JA,I1
               II=II+1
               I2=I2+1
  100    H(II)=H(II)+E2A(I2)-DE2A(I2)
C     CONTRIB D, CNDO.
         DO 110 I1=JC+1,JB
            II=(I1*(I1+1))/2
  110    H(II)=H(II)+E2A(1)-DE2A(1)
  120 CONTINUE
C
C   'SIZE' OF H IS NROW * NCOL
C
      RETURN
      END