      SUBROUTINE DERI22 (C,B,WORK,NORBS,FOC2,AB,MINEAR,FCI)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension C(NORBS,NORBS), B(*), WORK(NORBS,NORBS), FOC2(*),
     1AB(*), FCI(*)
************************************************************************
*  1) BUILD THE 2-ELECTRON FOCK MATRIX DEPENDING ON B AS FOLLOWS :
*     DP = C * SCALE*B * C' ...  DP DENSITY MATRIX 'DERIVATIVE',
*     FOC2 = 0.5 * TRACE ( DP * (2<J>-<K>) ) DONE IN FOCK2 & FOCK1.
*  2) HALF-TRANSFORM ONTO M.O. BASIS : DPT =  FOC2 * C
*     AND COMPUTE DIAGONAL BLOCKS ELEMENTS OF C' * FOC2, EXTRACTING
*     IN FCI ELEMENTS OVER C.I-ACTIVE M.O ONLY.
*  3) COMPUTE SUPERVECTOR AB = (DIAG + A) * B DEFINED BY THE MATRIX :
*     AB(I,J)= ( DIAG(I,J)*B(I,J)+DPT(I,J) )*SCALAR(I,J)  WITH I.GT.J,
*     DIAG(I,J)=(EIGS(I)-EIGS(J))/(O(J)-O(I)) >0, O OCCUPANCY NUMBERS,
*     EIGS EIGENVALUES OF FOCK OPERATOR WITH EIGENVECTORS C IN A.O.
*
*   INPUT
* C(NORBS,NORBS)   : M.O. EIGENVECTORS (COLUMNWISE).
* B(*)             : B SUPERVECTOR PACKED BY OFF-DIAGONAL BLOCKS, SCALED
* WORK(*)          : WORK AREA OF SIZE N*N.
* NORBS            : NUMBER OF M.O.S
* NELEC,NMOS       : LAST FROZEN CORE M.O. , C.I-ACTIVE BAND LENGTH.
*           IN common
* DIAG,SCALAR AS DEFINED IN 'DERI0'.
*   OUTPUT
* FOC2(*)       : 2-ELECTRON FOCK MATRIX, PACKED CANONICAL.
* AB(*)         : ANTISYMMETRIC MATRIX PACKED IN SUPERVECTOR FORM WITH
*                 THE CONSECUTIVE FOLLOWING BLOCKS:
*              1) OPEN-CLOSED  I.E. B(IJ)=B(I,J) WITH I OPEN & J CLOSED
*                 AND I RUNNING FASTER THAN J,
*              2) VIRTUAL-CLOSED SAME RULE OF ORDERING,
*              3) VIRTUAL-OPEN   SAME RULE OF ORDERING.
* FCI(*)        : FOCK DIAGONAL BLOCKS ELEMENTS OVER C.I-ACTIVE M.O.
*            FOC2 CAN BE equivalenceD WITH WORK IN THE callING SEQUENCE.
************************************************************************
C
C  NOTE: NORBS AND NORD ARE THE SAME ADDRESS.  THE NAME NORBD IS NOT
C        USED HERE.
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM)
     1               ,NLAST(NUMATM),NORBD,NELECS,NALPHA,NBETA
     2               ,NCLOSE,NOPEN,NDUMY,FRACT
     3       /WMATRX/ WJ(N2ELEC),WK(N2ELEC)
      common /DENSTY/ PDUMY(MPACK*2), DPA(MPACK)
      common /FOKMAT/ FDUMY(MPACK), SCALAR(MPACK)
      common /NVOMAT/ DIAG(MPACK/2)
      common /WORK1 / FDUMY2(15*MPACK), DP(7*MPACK)
      common /CIBITS/ NMOS,LAB,NELEC,NBO(3)
      dimension W(N2ELEC)
      equivalence (W,WJ)
C
      LINEAR=(NORBS*(NORBS+1))/2
C
C     DERIVATIVE OF THE DENSITY MATRIX IN DP (PACKED,CANONICAL).
C     ----------------------------------------------------------
C     DP = C * B * C' .
C
C     STEP 0 : UNSCALE VECTOR B.
      do 10 I=1,MINEAR
   10 B(I)=B(I)*SCALAR(I)
C
C     STEP 1 : WORK = C * B    .  DP TEMPORARY ARRAY.
      L=1
      if (NBO(2).NE.0 .AND. NBO(1).NE.0) then
C        OPEN-CLOSED
         call MXM(C(1,NBO(1)+1),NORBS,B(L),NBO(2),WORK,NBO(1))
C        CLOSED-OPEN
         call MXMT  (C,NORBS,B(L),NBO(1),WORK(1,NBO(1)+1),NBO(2))
         L=L+NBO(2)*NBO(1)
      end if
      if (NBO(3).NE.0 .AND. NBO(1).NE.0) then
C        VIRTUAL-CLOSED
         if (L.GT.1) then
            call MXM(C(1,NOPEN+1),NORBS,B(L),NBO(3),DP,NBO(1))
            do 20 I=1,NORBS*NBO(1)
   20       WORK(I,1)=WORK(I,1)+DP(I)
         else
            call MXM(C(1,NOPEN+1),NORBS,B(L),NBO(3),WORK,NBO(1))
         end if
C        CLOSED-VIRTUAL
         call MXMT(C,NORBS,B(L),NBO(1),WORK(1,NOPEN+1),NBO(3))
         L=L+NBO(3)*NBO(1)
      end if
      if (NBO(3).NE.0 .AND. NBO(2).NE.0) then
C        VIRTUAL-OPEN
         call MXM(C(1,NOPEN+1),NORBS,B(L),NBO(3),DP,NBO(2))
         J=NORBS*NBO(1)
         do 30 I=1,NORBS*NBO(2)
   30    WORK(J+I,1)=WORK(J+I,1)+DP(I)
C        OPEN-VIRTUAL
         call MXMT  (C(1,NBO(1)+1),NORBS,B(L),NBO(2),DP,NBO(3))
         J=NORBS*NOPEN
         do 40 I=1,NORBS*NBO(3)
   40    WORK(J+I,1)=WORK(J+I,1)+DP(I)
      end if
C
C     STEP 2 : DP= WORK * C'   WITH DP PACKED,CANONICAL.
      L=0
      do 50 I=1,NORBS
         do 50 J=1,I
            L=L+1
   50 DP(L)=SDOT(NORBS,WORK(I,1),NORBS,C(J,1),NORBS)
C
C     2-ELECTRON FOCK MATRIX BUILD WITH THE DENSITY MATRIX DERIVATIVE.
C     ----------------------------------------------------------------
C     RETURNED IN FOC2 (PACKED CANONICAL).
      do 60 I=1,LINEAR
         FOC2(I)=0.D0
   60 DPA(I)=0.5D0*DP(I)
      call FOCK2 (FOC2,DP,DPA,W,WJ,WK,NUMAT,NFIRST,NMIDLE,NLAST)
      call FOCK1 (FOC2,DP,DPA,DPA)
C
C     BUILD DP AND EXTRACT FCI.
C     --------------------------
C
C     DP(NORBS,Nend) = FOC2(NORBS,NORBS) * C(NORBS,NEND).
      Nend=MAX(NOPEN,NELEC+NMOS)
      L=1
      do 70 I=1,NOPEN
         call SUPDOT (DP(L),FOC2,C(1,I),NORBS,1)
   70 L=L+NORBS
C     EXTRACT FCI
      L=1
      Nend=0
      do 90 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 90
         do 80 I=N1,N2
            if (I.GT.NINIT) then
               call MXM (C(1,I),1,DP(NORBS*(NINIT-1)+1),NORBS,FCI(L),I-N
     1INIT)
               L=L+I-NINIT
            end if
   80    continue
   90 continue
      do 100 I=NELEC+1,NELEC+NMOS
         FCI(L)=-DOT(C(1,I),DP(NORBS*(I-1)+1),NORBS)
  100 L=L+1
C
C     NEW SUPERVECTOR AB = (DIAG + C'* FOC2 * C) * B , SCALED.
C     --------------------------------------------------------
C
C     PART 1 : AB(I,J) = (C' * DP) (I,J) DONE BY BLOCKS.
      L=1
      if (NBO(2).NE.0 .AND. NBO(1).NE.0) then
         call MTXM (C(1,NBO(1)+1),NBO(2),DP,NORBS,AB(L),NBO(1))
         L=L+NBO(2)*NBO(1)
      end if
      if (NBO(3).NE.0 .AND. NBO(1).NE.0) then
         call MTXM (C(1,NOPEN+1),NBO(3),DP,NORBS,AB(L),NBO(1))
         L=L+NBO(3)*NBO(1)
      end if
      if (NBO(3).NE.0 .AND. NBO(2).NE.0)
     1    call MTXM(C(1,NOPEN+1),NBO(3),DP(NORBS*NBO(1)+1),
     2                                   NORBS,AB(L),NBO(2))
C
C     PART 2 : AB = SCALE * (D * B + AB) AND RESCALE BASIS VECTOR B.
      do 110 I=1,MINEAR
         AB(I)=(DIAG(I)*B(I)+AB(I))*SCALAR(I)
  110 B(I)=B(I)/SCALAR(I)
      return
      end
