      SUBROUTINE MECIP(COEFFS,NORBS,DELTAP, DELTA)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension COEFFS(NORBS,NORBS), DELTAP(NMOS,NMOS),
     1 DELTA(NORBS,NMOS)
************************************************************************
*
*   MECIP WILL CORRECT THE TOTAL DENSITY MATRIX FOR THE EFFECT OF THE
*   C.I.
*              ON INPUT
*
*  COEFFS       : ALL M.O.'S (NORBS M.O.S)
*  NORBS        : NUMBER OF MOLECULAR ORBITALS = NUMBER OF A.O.'S
*  P            : TOTAL DENSITY MATRIX
*  NMOS         : NUMBER OF M.O.'S IN ACTIVE SPACE
*  VECTCI       : STATE VECTOR OF LENGTH LAB
*  MICROA(I,J)  : ALPHA OCCUPANCY OF M.O. 'I' IN MICROSTATE 'J'
*  MICROB(I,J)  : BETA  OCCUPANCY OF M.O. 'I' IN MICROSTATE 'J'
*
*  NOTE: THIS IS A MODIFICATION OF CODE ORIGINALLY WRITTEN BY
*        PROF. DANIEL LIOTARD
************************************************************************
      common /CIBITS/ NMOS,LAB,NELEC, NBO(3)
      common /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
      common /NALMAT/ NALPHA(NMECI**2)
      common /BASEOC/ OCCA(NMECI)
      common /CIVECT/ VECTCI(NMECI**2),CONF(NMECI**4)
      common /MICROS/ MICROA(NMECI,4*NMECI**2), MICROB(NMECI,4*NMECI**2)
C     INITIALIZE WITH THE OPPOSITE OF THE 'SCF' DENSITY.
      do 10 I=1,NMOS
         DELTAP(I,I)=-OCCA(I)*2.D0
         do 10 J=1,I-1
   10 DELTAP(I,J)=0.D0
C
C     ADD THE C.I. CORRECTION
      do 120 ID=1,LAB
         do 120 JD=1,ID
C     CHECK SPIN AGREEMENT
            if (NALPHA(ID).NE.NALPHA(JD)) go to 120
            IX=0
            IY=0
            do 20 J=1,NMOS
               IX=IX+ABS(MICROA(J,ID)-MICROA(J,JD))
   20       IY=IY+ABS(MICROB(J,ID)-MICROB(J,JD))
C     CHECK NUMBER OF DIFFERING M.O.
            if (IX+IY.GT.2) go to 120
            if (IX.EQ.2) then
C        DETERMINANTS ID AND JD DIFFER BY M.O I IN ID AND M.O J IN JD:
               do 30 I=1,NMOS
   30          if (MICROA(I,ID).NE.MICROA(I,JD)) go to 40
   40          IJ=MICROB(I,ID)
               do 50 J=I+1,NMOS
                  if (MICROA(J,ID).NE.MICROA(J,JD)) go to 60
   50          IJ=IJ+MICROA(J,ID)+MICROB(J,ID)
C        IJ GIVES THE SIGN OF THE PERMUTATION
   60          DELTAP(J,I)=DELTAP(J,I)+VECTCI(ID)*VECTCI(JD)*FLOAT(1-2*M
     1OD(IJ,2))
            else if (IY.EQ.2) then
C        DETERMINANTS ID AND JD DIFFER BY M.O J IN ID AND M.O I IN JD:
               do 70 I=1,NMOS
   70          if (MICROB(I,ID).NE.MICROB(I,JD)) go to 80
   80          IJ=0
               do 90 J=I+1,NMOS
                  if (MICROB(J,ID).NE.MICROB(J,JD)) go to 100
   90          IJ=IJ+MICROA(J,ID)+MICROB(J,ID)
  100          IJ=IJ+MICROA(J,ID)
               DELTAP(J,I)=DELTAP(J,I)+VECTCI(ID)*VECTCI(JD)*FLOAT(1-2*M
     1OD(IJ,2))
            else
C        DETERMINANTS ID AND JD ARE IDENTICAL:
               do 110 I=1,NMOS
  110          DELTAP(I,I)=DELTAP(I,I)+(MICROA(I,ID)+MICROB(I,ID))*
     1                                                VECTCI(ID)**2
            end if
  120 continue
C
C     BACK TRANSFORM INTO A.O. BASIS.
C     -------------------------------
C     P(C.I.) = P(SCF) + C * DELTAP * C'
      do 130 I=1,NMOS
CDIR$ IVDEP
         do 130 J=1,I-1
  130 DELTAP(J,I)=DELTAP(I,J)
C     STEP 1: DELTAP = C * DELTAP
      call MXM (COEFFS(1,NELEC+1),NORBS,DELTAP,NMOS,DELTA,NMOS)
C     STEP 2: P = P + DELTAP * C'
      IJ=0
      do 150 I=1,NORBS
         do 150 J=1,I
            IJ=IJ+1
            SUM=0.D0
            do 140 K=1,NMOS
  140       SUM=SUM+DELTA(I,K)*COEFFS(J,NELEC+K)
  150 P(IJ)=P(IJ)+SUM
C     NOTE FROM D.L.: AT THIS POINT THE 'NATURAL ORBITALS' OF THIS STATE
C     CAN BE OBTAINED STRAIGHTWAY AS EIGENVECTORS OF THE DENSITY MATRIX.
      return
      end
