      SUBROUTINE DERI2(C,E,NORBS, MINEAR, F, FD, FCI, NINEAR,
     1NVAR,WORK,B,NW2,GRAD,AB,NW3,FB,THROLD)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension C(NORBS,NORBS),E(NORBS),WORK(NORBS,NORBS)
     1         ,F(MINEAR,MPACK/MINEAR)
     2         ,FD(NINEAR,(5*MAXPAR+MPACK)/NINEAR)
     3         ,FCI(NINEAR,MAXORB/NINEAR)
     4         ,B(MINEAR,(6*MPACK-5*MAXPAR)/MINEAR)
     5         ,AB(MINEAR,6*MPACK/MINEAR),FB(NVAR,6*MPACK/NVAR)
     6         ,GRAD(NVAR)
*********************************************************************
*
*     DERI2 COMPUTE THE RELAXATION PART OF THE DERIVATIVES OF THE
*     NON-VARIATIONALLY OPTIMIZED ENERGY WITH RESPECT TO TWO
*     COORDINATES AT A TIME. THIS IS DONE IN THREE STEPS.
*
*     THE M.O DERIVATIVES ARE SOLUTION {X} OF A LINEAR SYSTEM
*                        (D-A) * X = F
*     WHERE D IS A DIAGONAL SUPER-MATRIX OF FOCK EIGENVALUE DIFFERENCES
*     AND A IS A SUPER-MATRIX OF 2-ELECTRONS INTEGRALS IN M.O BASIS.
*     SUCH A SYSTEM IS TOO LARGE TO BE INVERTED DIRECTLY THUS ONE MUST
*     USES A RELAXATION METHOD TO GET A REASONABLE ESTIMATE OF {X}.
*     THIS REQUIRES A BASIS SET {B} TO BE GENERATED ITERATIVELY, AFTER
*     WHICH WE SOLVE BY DIRECT INVERSION THE LINEAR SYSTEM PROJECTED
*     IN THIS BASIS {B}. IT WORKS QUICKLY BUT DOES REQUIRE A LARGE
*     CORE MEMORY.
*
*     USE A FORMALISM WITH FOCK OPERATOR THUS AVOIDING THE EXPLICIT
*     COMPUTATION (AND STORAGE) OF THE SUPER-MATRIX A.
*     THE SEMIEMPIRICAL METHODS do NOT INVOLVE LARGE C.I CALCULATIONS.
*     THEREFORE FOR EACH GRADIENT COMPONENT WE BUILD THE C.I MATRIX
*     DERIVATIVE FROM THE M.O. INTEGRALS <IJ|KL> AND FOCK EIGENVALUES
*     DERIVATIVES, THUS PROVIDING THE RELAXATION CONTRIBUTION TO THE
*     GRADIENT WITHOUT COMPUTATION AND STORAGE OF THE 2ND ORDER DENSITY
*     MATRIX.
*
*   STEP 1)
*     USE THE PREVIOUS B AND THE NEW F VECTORS TO BUILD AN INITIAL
*     BASIS SET B.
*   STEP 2)
*     BECAUSE THE ELECTRONIC HESSIAN (D-A) IS THE SAME FOR EACH
*     DERIVATIVE, WE ONLY NEED TO ENLARGE ITERATIVELY THE ORTHONORMAL
*     BASIS SET {B} USED TO INVERT THE PROJECTED HESSIAN.
*     (DERIVED FROM THE LARGEST RESIDUAL VECTOR ).
*     THIS SECTION IS CARRIED OUT IN THE DIAGONAL METRIC 'SCALAR'.
*   STEP 3) ... LOOP ON THE GEOMETRIC VARIABLE :
* 3.1 FOR EACH GEOMETRIC VARIABLE, GET THE M.O DERIVATIVES IN A.O.
* 3.2 COMPUTE THE FOCK EIGENVALUES AND 2-ELECTRON INTEGRAL RELAXATION.
* 3.3 BUILD THE ELECTRONIC RELAXATION CONTRIBUTION TO THE C.I MATRIX
*     AND GET THE ASSOCIATED EIGENSTATE DERIVATIVE WITH RESPECT TO
*     THE GEOMETRIC VARIABLE.
*
*   INPUT
*     C(NORBS,NORBS) : M.O. COEFFICIENTS, IN COLUMN.
*     E(NORBS)       : EIGENVALUES OF THE FOCK MATRIX.
*     MINEAR         : NUMBER OF NON REDUNDANT ROTATION OF THE M.O.
*     F(MINEAR,NVAR) : NON-RELAXED FOCK MATRICES DERIVATIVES
*                    IN M.O BASIS, OFF-DIAGONAL BLOCKS.
*     FD(NINEAR,NVAR): IDEM, DIAGONAL BLOCKS, C.I-ACTIVE ONLY.
*     WORK           : WORK ARRAY OF SIZE N*N.
*     B(MINEAR,NBSIZE) : INITIAL ORTHONORMALIZED BASIS SET {B}.
*     GRAD(NVAR)     : GRADIENT VECTOR BEFORE RELAXATION CORRECTION.
*     AB(MINEAR,*): STORAGE FOR THE (D-A) * B VECTORS.
*     FB(NVAR,*)  : STORAGE FOR THE MATRIX PRODUCT F' * B.
*   OUTPUT
*     GRAD   : DERIVATIVE OF THE HEAT OF FORMATION WITH RESPECT TO
*              THE NVAR OPTIMIZED VARIABLES.
*
************************************************************************
      common /FOKMAT/ FDUMY(MPACK), SCALAR(MPACK)
      common /NVOMAT/ DIAG(MPACK/2)
     1       /WORK3 / DIJKL(MPACK*4)
     2       /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
     3       /CIVECT/ VECTCI(NMECI**2),BABINV(NMECI**3),
     4BCOEF(NMECI**4-NMECI**3)
     5       /KEYWRD/ KEYWRD
      common /CIBITS/ NMOS,LAB,NELEC,NBO(3)
      common /WORK2 / BAB(60,60),DUMY(NMECI**4+2*NMECI**3+NMECI**2-3600)
      common /NUMCAL/ NUMCAL
      dimension LCONV(60)
      logical FAIL, LCONV, DEBUG, LBAB
      character KEYWRD*241
      data ICALCN/0/
C
C     * * * STEP 1 * * *
C     BUILD UP THE INITIAL ORTHONORMALIZED BASIS.
C
      if (ICALCN.NE.NUMCAL) then
         DEBUG=INDEX(KEYWRD,' DERI2').NE.0
         ICALCN=NUMCAL
         LINEAR=NORBS*(NORBS+1)/2
         NOPEN =NBO(1)+NBO(2)
         MAXITE=MIN(60,INT(SQRT(NMECI**3.D0)),MPACK*2/NVAR)
         MAXITE=MIN(MAXITE,MIN(NW2,NW3)/MAX(MINEAR,NINEAR))
         NFIRST=MIN(NVAR,1+MAXITE/4)
      end if
      FAIL=.FALSE.
      NBSIZE=0
      TIME1=SECOND()
C
C        NORMAL CASE. USE F ONLY.
C
      call DERI21 (F,NVAR,MINEAR,NFIRST,WORK
     1               ,WORK(NVAR*NVAR+1,1),B,NLAST)
      LBAB=.FALSE.
      NFIRST=NBSIZE+1
      NLAST=NBSIZE+NLAST
      do 10 I=1,NVAR
   10 LCONV(I)=.FALSE.
C
C     * * * STEP 2 * * *
C     RELAXATION METHOD WITH OPTIMUM INCREASE OF THE BASIS SET.
C     ---------------------------------------------------------
C
C     UPDATE AB ,FCI AND BAB. (BAB IS SYMMETRIC)
   20 do 30 J=NFIRST,NLAST
         call DERI22(C,B(1,J),WORK,NORBS,WORK,AB(1,J),MINEAR,
     1            FCI(1,J))
         call MXM(AB(1,J),1,B,MINEAR,BAB(1,J),NLAST)
         do 30 I=1,NFIRST-1
   30 BAB(J,I)=BAB(I,J)
C     INVERT BAB, STORE IN BABINV.
   40 L=0
      do 50 J=1,NLAST
         do 50 I=1,NLAST
            L=L+1
   50 BABINV(L)=BAB(I,J)
      call OSINV (BABINV,NLAST,DETER)
      if (DETER.EQ.0) then
      if (NLAST.NE.1) then
         write(MIOB,'('' THE BAB MATRIX OF ORDER'',I3,
     1   '' IS SINGULAR IN DERI2''/
     2   '' THE RELAXATION IS STOPPED AT THIS POINT.'')') NLAST
         end if
         LBAB=.TRUE.
         NLAST=NLAST-1
         go to 40
      end if
      if (.NOT.LBAB) then
C        UPDATE F * B'
         call MTXM (F,NVAR,B(1,NFIRST),MINEAR,FB(1,NFIRST),NLAST-NFIRST+
     11)
      end if
C     NEW SOLUTIONS IN BASIS B , STORED IN BCOEF(NVAR,*).
C     BCOEF = BABINV * FB'
      if (NLAST.NE.0) call MXMT (BABINV,NLAST,FB,NLAST,BCOEF,NVAR)
      if (LBAB) go to 90
C
C     SELECT THE NEXT BASIS VECTOR AS THE LARGEST RESIDUAL VECTOR.
C     AND TEST FOR CONVERGENCE ON THE LARGEST RESIDUE.
      NRES=0
      TEST2=0.D0
      do 70 IVAR=1,NVAR
         if (LCONV(IVAR)) go to 70
C     GET ONE NOT-CONVERGED RESIDUAL VECTOR (# IVAR),
C     STORED IN WORK.
         call MXM  (AB,MINEAR,BCOEF(NLAST*(IVAR-1)+1),NLAST,WORK,1)
         TEST=0.D0
         do 60 I=1,MINEAR
            WORK(I,1)=F(I,IVAR)-WORK(I,1)
   60    TEST=MAX(ABS(WORK(I,1)),TEST)
         if (DEBUG) write(MIOB,*) ' TEST:', TEST
         TEST2=MAX(TEST2,TEST)
         if (TEST.LE.THROLD) then
            LCONV(IVAR)=.TRUE.
            if (NVAR.EQ.1) go to 90
            go to 70
         else if (NLAST+NRES.EQ.MAXITE-1) then
C        RUNNING OUT OF STORAGE
            if (TEST.LE.MAX(0.01D0,THROLD*2)) then
               LCONV(IVAR)=.TRUE.
               go to 70
            end if
         else if (NLAST+NRES.EQ.MAXITE) then
C
C   COMPLETELY OUT OF STORAGE
C
            FAIL=NRES.EQ.0
            go to 80
         else
C        STORE THE FOLLOWING RESIDUE IN AB(continueD).
            NRES=NRES+1
            call SCOPY (MINEAR,WORK,1,AB(1,NLAST+NRES),1)
         end if
   70 continue
   80 if (NRES.EQ.0) go to 90
C     FIND OPTIMUM FOLLOWING SUBSET, ADD TO B AND LOOP.
      NFIRST=NLAST+1
      call DERI21(AB(1,NFIRST),NRES,MINEAR,NRES,WORK
     1           ,WORK(NRES*NRES+1,1),B(1,NFIRST),NADD)
      NLAST=NLAST+NADD
      go to 20
C
C     CONVERGENCE ACHIEVED OR HALTED.
C     -------------------------------
C
   90 NBSZE=NBSIZE
      if (DEBUG.OR.LBAB) then
         write(MIOB,'('' RELAXATION endED IN DERI2 AFTER'',I3,
     1   '' CYCLES''/'' REQUIRED CONVERGENCE THRESHOLD ON RESIDUALS =''
     2   ,F12.9/'' HIGHEST RESIDUAL ON'',I3,'' GRADIENT COMPONENTS = ''
     3   ,F12.9)') NLAST-NBSZE, THROLD, NVAR, TEST2
      if (NLAST-NBSZE.EQ.0) then
      write(MIOB,'(A)')
     +' ANALYTIC C.I. DERIVATIVES do NOT WORK FOR THIS SYSTEM'
      write(MIOB,'(A)') ' ADD KEYWORD ''NOANCI'' AND RESUBMIT'
      call filclo  ! J/M
      stop ' in DERI2 - analytic CI derivs. not applicable'
      end if
         TIME2=SECOND()
         write(MIOB,
     1          '('' ELAPSED TIME IN RELAXATION'',F15.3,'' SECOND'') ')
     2   TIME2-TIME1
      end if
      if (FAIL) then
         write(MIOB,'(A)')
     1              ' ANALYTICAL DERIVATIVES TOO INACCURATE FOR THIS'
         write(MIOB,'(A)')
     1              ' WORK.  JOB STOPPED HERE.  SEE MANUAL FOR IDEAS'
         call filclo  ! J/M
         stop ' in DERI2 - analytic derivs. not accurate enough'
      else
         NBSIZE=0
C        UNSCALED SOLUTION SUPERVECTORS, STORED IN F.
         if (NLAST.NE.0) call MXM (B,MINEAR,BCOEF,NLAST,F,NVAR)
         do 100 J=1,NVAR
            do 100 I=1,MINEAR
  100    F(I,J)=F(I,J)*SCALAR(I)
C        FOCK MATRIX DIAGONAL BLOCKS OVER C.I-ACTIVE M.O.
C        STORED IN FB.
         if (NLAST.NE.0) call MXM (FCI,NINEAR,BCOEF,NLAST,FB,NVAR)
      end if
C
C     * * * STEP 3 * * *
C     FINAL LOOP (390) ON THE GEOMETRIC VARIABLES.
C     --------------------------------------------
C
      do 130 IVAR=1,NVAR
C
C     C.I-ACTIVE M.O DERIVATIVES INTO THE M.O BASIS,
C         RETURNED IN AB (N,NELEC+1,...,NELEC+NMOS).
C     C.I-ACTIVE EIGENVALUES DERIVATIVES,
C         RETURNED IN BCOEF(NELEC+1,...,NELEC+NMOS).
         call DERI23 (F(1,IVAR),MINEAR,FD(1,IVAR),NINEAR,E
     1            ,FB(NINEAR*(IVAR-1)+1,1),AB,BCOEF,NORBS)
C
C     DERIVATIVES OF THE 2-ELECTRONS INTEGRALS OVER C.I-ACTIVE M.O.
C     STORED IN /XYIJKL/.
         call DIJKL2 (AB(NORBS*NELEC+1,1),NORBS,NMOS,DIJKL,XY,NMECI)
         if (DEBUG) then
            write(MIOB,'('' * * * GRADIENT COMPONENT NUMBER'',I4)') IVAR
            if (INDEX(KEYWRD,'DEBU').NE.0) then
               write(MIOB,
     1                   '('' C.I-ACTIVE M.O. DERIVATIVES IN M.O BASIS'',
     2                     '', IN ROW.'')')
               L=NORBS*NELEC+1
               do 110 I=NELEC+1,NELEC+NMOS
                  write(MIOB,'(8F10.4)') (AB(K,1),K=L,L+NORBS-1)
  110          L=L+NORBS
            end if
            write(MIOB,
     1       '('' C.I-ACTIVE FOCK EIGENVALUES RELAXATION (E.V.)'')')
            write(MIOB,'(8F10.4)') (BCOEF(I),I=NELEC+1,NELEC+NMOS)
            write(MIOB,'('' 2-ELECTRON INTEGRALS RELAXATION (E.V.)''/
     1''    I    J    K    L       d<I(1)J(1)|K(2)L(2)> RELAXATION ONLY'
     2')
     3')
            do 120 I=1,NMOS
               do 120 J=1,I
                  do 120 K=1,I
                     LL=K
                     if (K.EQ.I) LL=J
                     do 120 L=1,LL
  120       write(MIOB,'(4I5,F20.10)')
     1              NELEC+I,NELEC+J,NELEC+K,NELEC+L,XY(I,J,K,L)
         end if
C
C     BUILD THE C.I MATRIX DERIVATIVE, STORED IN AB.
         call MECID (BCOEF,GSE,WORK(LAB+1,1),WORK)
         call MECIH (WORK,AB,NMOS,LAB)
C     RELAXATION CORRECTION TO THE C.I ENERGY DERIVATIVE.
         call SUPDOT (WORK,AB,VECTCI,LAB,1)
         GRAD(IVAR)=GRAD(IVAR)+DOT(VECTCI,WORK,LAB)*23.061D0
         if (DEBUG) then
            write(MIOB,'('' RELAXATION OF THE GRADIENT COMPONENT'',
     1       F10.4,'' KCAL/MOLE'')')  DOT(VECTCI,WORK,LAB)*23.061D0
         end if
C
C     THE end .
  130 continue
      if (DEBUG)
     1  write(MIOB,'('' ELAPSED TIME IN C.I-ENERGY RELAXATION'',F15.3,
     2               '' SECOND'')') SECOND()-TIME2
      return
      end
