      SUBROUTINE IJKL(CP, CF, NELEC, NMOS, DIJKL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION CP(NORBS,NMOS), CF(NORBS,NORBS),
     1DIJKL(NORBS,NMOS,(NMOS*(NMOS+1))/2)
************************************************************************
*
*   IJKL FILLS THE XY ARRAY.  XY HOLDS THE TWO-ELECTRON INTEGRALS OVER
*        MOLECULAR ORBITALS IN THE ACTIVE SPACE.
*        XY(I,J,K,L) = <IJ|1/R(1,2)|KL>
*
*           ON INPUT
*
* CP     = M.O.'S OVER C.I. ACTIVE SPACE (NORMALLY 1 TO 5 M.O.S)
* CF     = ALL M.O.'S, INCLUDING THOSE IN CP
* NORBS  = NUMBER OF ATOMIC ORBITALS
* NELEC  = NUMBER OF OCCUPIED M.O.S NOT INVOLVED IN THE C.I.
* NMOS   = NUMBER OF M.O.S INVOLVED IN THE C.I. (NORMALLY 1 TO 5 M.O.S)
*          ALSO CALLED THE ACTIVE SPACE OF THE C.I.
*
*  NOTE: THIS SUBROUTINE IS UNUSUAL IN THAT ONE FUNCTION IS TO
*        FILL THE ARRAY XY WHICH IS NOT PASSED AS AN ARGUMENT,
*        INSTEAD IT IS PASSED VIA COMMON BLOCK XYIJKL.
*
************************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /VECTOR/ DUMY1(MORB2+MAXORB), CIJ(10*MAXORB),
     1WCIJ(10*MAXORB), DUMY(MORB2-19*MAXORB)
      COMMON /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
      DIMENSION CKL(10*MAXORB)
      DIMENSION NB(0:8)
      EQUIVALENCE (CKL, CIJ)
      DATA NB /1,0,0,10,0,0,0,0,45/
C
C  CALCULATE TWO-ELECTRON INTEGRALS FOR THE SET DIJKL(K,L,IJ)
C  THE INDEX K RUNS OVER ALL M.O.'S, L OVER ACTIVE-SPACE M.O.'S,
C  AND IJ OVER LOWER-HALF TRIANGLE OF ACTIVE-SPACE M.O.'S, J FASTER THAN
C  I.
C  ALL ACTIVE-SPACE INTERACTIONS ARE COPIED INTO THE ARRAY XY
C
      IJ=0
      DO 70 I=1,NMOS
         DO 70 J=1,I
            IJ=IJ+1
            IPQ=0
            DO 10 II=1,NUMAT
               DO 10 IP=NFIRST(II),NLAST(II)
                  DO 10 IQ=NFIRST(II),IP
                     IPQ=IPQ+1
   10       CIJ(IPQ)=CP(IP,I)*CP(IQ,J)+CP(IP,J)*CP(IQ,I)
C
C  CIJ HOLDS THE DENSITY DISTRIBUTION PSI(I)*PSI(J) OVER ATOMIC BASES
C  I AND J ARE M.O. INDICES WITHIN THE ACTIVE SPACE.  CIJ(M,N) IS FOR
C  THE ATOMIC BASES M AND N FOR M.O.'S I AND J.
C
            CALL PARTXY (CIJ,WCIJ)
C
C WCIJ HOLDS THE KET PART OF THE INTEGRAL <K,L|1/R(1,2)|I,J>
C THAT IS, |1/R(1,2)|I,J>.  WCIJ(M,N) IS FOR THE ATOMIC BASES M AND N
C FOR M.O.'S K AND L.
C
   20       CONTINUE
            DO 60 K=1,NORBS
               DO 50 L=1,NMOS
C
C  ABOUT TO CALCULATE <I,J|1/R(1,2)|K,L>
C
                  IPQ=0
                  DO 30 II=1,NUMAT
                     DO 30 IP=NFIRST(II),NLAST(II)
                        DO 30 IQ=NFIRST(II),IP
                           IPQ=IPQ+1
   30             CKL(IPQ)=CF(IP,K)*CP(IQ,L)+CP(IP,L)*CF(IQ,K)
C
C CKL HOLDS THE DENSITY DISTRIBUTION PSI(K)*PSI(L) OVER ATOMIC BASES.
C K IS THE INDEX OF A M.O.; L IS AN INDEX OF A M.O. IN THE ACTIVE SPACE.
C
                  SUM=0.D0
                  DO 40 II=1,IPQ
   40             SUM=SUM+CKL(II)*WCIJ(II)
C
C  SUM IS THE INTEGRAL <I,J|1/R(1,2)|K,L>
C
                  DIJKL(K,L,IJ)=SUM
   50          CONTINUE
   60       CONTINUE
   70 CONTINUE
C
C  NOW SPREAD THE INTEGRALS OVER THE XY ARRAY.  XY IS ENTIRELY
C  IN ACTIVE SPACE
C
      DO 80 K=1,NMOS
         KK=NELEC+K
C
C  K IS A M.O. INDEX IN ACTIVE SPACE
C KK IS A M.O. INDEX
C
         DO 80 L=1,NMOS
            IJ=0
            DO 80 I=1,NMOS
               DO 80 J=1,I
                  IJ=IJ+1
                  SUM=DIJKL(KK,L,IJ)
                  XY(I,J,K,L)=SUM
                  XY(I,J,L,K)=SUM
                  XY(J,I,K,L)=SUM
                  XY(J,I,L,K)=SUM
                  XY(K,L,I,J)=SUM
                  XY(K,L,J,I)=SUM
                  XY(L,K,I,J)=SUM
   80 CONTINUE
      RETURN
      END
