      SUBROUTINE FOCK2(F, PTOT, P, W, WJ, WK, NUMAT, NFIRST,
     1NMIDLE, NLAST)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension F(*), PTOT(*), WJ(*), WK(*), NFIRST(*), NMIDLE(*),
     1          NLAST(*), P(*), W(*)
      double precision WJ,WK
C***********************************************************************
C
C FOCK2 FORMS THE TWO-ELECTRON TWO-CENTER REPULSION PART OF THE FOCK
C MATRIX
C ON INPUT  PTOT = TOTAL DENSITY MATRIX.
C           P    = ALPHA OR BETA DENSITY MATRIX.
C           W    = TWO-ELECTRON INTEGRAL MATRIX.
C
C  ON OUTPUT F   = PARTIAL FOCK MATRIX
C***********************************************************************
      common /EULER / TVEC(3,3), ID
      common /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
      common /NUMCAL/ NUMCAL
      common /KEYWRD/ KEYWRD
      common /WORK4 / PTOT2
      dimension IFACT(MAXORB),
     1I1FACT(MAXORB), JINDEX(256), KINDEX(256), IJPERM(10), LLPERM(10),
     2PK(16), PJA(16), PJB(16), MMPERM(10),
     3PTOT2(MAXORB,MAXORB), JJNDEX(256)
      character*241 KEYWRD
      logical LID
      save IFACT,I1FACT, IONE, LID
      data ICALCN/0/
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
C
C   SET UP ARRAY OF LOWER HALF TRIANGLE INDICES (PASCAL'S TRIANGLE)
C
         do 10 I=1,MAXORB
            IFACT(I)=(I*(I-1))/2
   10    I1FACT(I)=IFACT(I)+I
C
C   SET UP GATHER-SCATTER TYPE ARRAYS FOR USE WITH TWO-ELECTRON
C   INTEGRALS.  JINDEX ARE THE INDICES OF THE J-INTEGRALS FOR ATOM I
C   INTEGRALS.  JJNDEX ARE THE INDICES OF THE J-INTEGRALS FOR ATOM J
C               KINDEX ARE THE INDICES OF THE K-INTEGRALS
C
         M=0
         do 20 I=1,4
            do 20 J=1,4
               IJ=MIN(I,J)
               JI=I+J-IJ
               do 20 K=1,4
                  IK=MIN(I,K)
                  KI=I+K-IK
                  do 20 L=1,4
                     M=M+1
                     KL=MIN(K,L)
                     LK=K+L-KL
                     JL=MIN(J,L)
                     LJ=J+L-JL
                     KINDEX(M)= IFACT(LJ) +JL + 10*( IFACT(KI) +IK) -10
   20    JINDEX(M)=(IFACT(JI) + IJ)*10 + IFACT(LK) + KL - 10
         L=0
         do 30 I=1,4
            I1=(I-1)*4
            do 30 J=1,I
               I1=I1+1
               L=L+1
               IJPERM(L)=I1
               MMPERM(L)=IJPERM(L)-16
               LLPERM(L)=(I1-1)*16
   30    continue
         L=0
         do 40 I=1,10
            M=MMPERM(I)
            L=LLPERM(I)
            do 40 K=1,16
               L=L+1
               M=M+16
   40    JJNDEX(L)=JINDEX(M)
         LID=(ID.EQ.0)
         IONE=1
         if (ID.NE.0) IONE=0
C
C      end OF INITIALIZATION
C
      end if
      if (ITYPE.EQ.4) go to 200
C
C     START OF MNDO, AM1, OR PM3 OPTION
C
      KK=0
      L=0
      do 60 I=1,NUMAT
         IA=NFIRST(I)
         IB=NLAST(I)
         M=0
         do 50 J=IA,IB
            do 50 K=IA,IB
               M=M+1
               JK=MIN(J,K)
               KJ=K+J-JK
               JK=JK+(KJ*(KJ-1))/2
               PTOT2(I,M)=PTOT(JK)
   50    continue
   60 continue
      do 190 II=1,NUMAT
         IA=NFIRST(II)
         IB=NLAST(II)
C
C  IF NUMAT=2 then WE ARE IN A DERIVATIVE OR IN A MOLECULE CALCULATION
C
         if (NUMAT.NE.2) then
            IMINUS=II-IONE
         else
            IMINUS=II-1
         end if
         do 180 JJ=1,IMINUS
            JA=NFIRST(JJ)
            JB=NLAST(JJ)
            JC=NMIDLE(JJ)
            if (LID) then
               DREP=W(KK+1)
               if (IB-IA.GE.3.AND.JB-JA.GE.3) then
C
C                         HEAVY-ATOM  - HEAVY-ATOM
C
C   EXTRACT COULOMB TERMS
C
                  do 70 I=1,16
                     PJA(I)=PTOT2(II,I)
   70             PJB(I)=PTOT2(JJ,I)
C
C  COULOMB TERMS
C
                  call JAB(IA,JA,LLPERM,JINDEX, JJNDEX, PJA,PJB,W(KK+1),
     1F)
C
C  EXCHANGE TERMS
C
C
C  EXTRACT INTERSECTION OF ATOMS II AND JJ IN THE SPIN DENSITY MATRIX
C
                  L=0
                  do 80 I=IA,IB
                     I1=IFACT(I)+JA
                     do 80 J=I1,I1+3
                        L=L+1
   80             PK(L)=P(J)
                  call KAB(IA,JA, PK, W(KK+1), KINDEX, F)
                  KK=KK+100
               else if (IB-IA.GE.3.AND.JA.EQ.JB) then
C
C                         LIGHT-ATOM  - HEAVY-ATOM
C
C
C   COULOMB TERMS
C
                  SUMDIA=0.D0
                  SUMOFF=0.D0
                  LL=I1FACT(JA)
                  K=0
                  do 100 I=0,3
                     J1=IFACT(IA+I)+IA-1
                     do 90 J=0,I-1
                        K=K+1
                        J1=J1+1
                        F(J1)=F(J1)+PTOT(LL)*W(KK+K)
   90                SUMOFF=SUMOFF+PTOT(J1)*W(KK+K)
                     J1=J1+1
                     K=K+1
                     F(J1)=F(J1)+PTOT(LL)*W(KK+K)
  100             SUMDIA=SUMDIA+PTOT(J1)*W(KK+K)
                  F(LL)=F(LL)+SUMOFF*2.D0+SUMDIA
C
C  EXCHANGE TERMS
C
C
C  EXTRACT INTERSECTION OF ATOMS II AND JJ IN THE SPIN DENSITY MATRIX
C
                  K=0
                  do 120 I=IA,IB
                     I1=IFACT(I)+JA
                     SUM=0.D0
                     do 110 J=IA,IB
                        K=K+1
                        J1=IFACT(J)+JA
  110                SUM=SUM+P(J1)*W(KK+JINDEX(K))
  120             F(I1)=F(I1)-SUM
                  KK=KK+10
               else if (JB-JA.GE.3.AND.IA.EQ.IB) then
C
C                         HEAVY-ATOM - LIGHT-ATOM
C
C
C   COULOMB TERMS
C
                  SUMDIA=0.D0
                  SUMOFF=0.D0
                  LL=I1FACT(IA)
                  K=0
                  do 140 I=0,3
                     J1=IFACT(JA+I)+JA-1
                     do 130 J=0,I-1
                        K=K+1
                        J1=J1+1
                        F(J1)=F(J1)+PTOT(LL)*W(KK+K)
  130                SUMOFF=SUMOFF+PTOT(J1)*W(KK+K)
                     J1=J1+1
                     K=K+1
                     F(J1)=F(J1)+PTOT(LL)*W(KK+K)
  140             SUMDIA=SUMDIA+PTOT(J1)*W(KK+K)
                  F(LL)=F(LL)+SUMOFF*2.D0+SUMDIA
C
C  EXCHANGE TERMS
C
C
C  EXTRACT INTERSECTION OF ATOMS II AND JJ IN THE SPIN DENSITY MATRIX
C
                  K=IFACT(IA)+JA
                  J=0
                  do 160 I=K,K+3
                     SUM=0.D0
                     do 150 L=K,K+3
                        J=J+1
  150                SUM=SUM+P(L)*W(KK+JINDEX(J))
  160             F(I)=F(I)-SUM
                  KK=KK+10
               else if (JB.EQ.JA.AND.IA.EQ.IB) then
C
C                         LIGHT-ATOM - LIGHT-ATOM
C
                  I1=I1FACT(IA)
                  J1=I1FACT(JA)
                  IJ=I1+JA-IA
                  F(I1)=F(I1)+PTOT(J1)*W(KK+1)
                  F(J1)=F(J1)+PTOT(I1)*W(KK+1)
                  F(IJ)=F(IJ)-P   (IJ)*W(KK+1)
                  KK=KK+1
               end if
            else
               DREP=WJ(KK+1)
               do 170 I=IA,IB
                  KA=IFACT(I)
                  do 170 J=IA,I
                     KB=IFACT(J)
                     IJ=KA+J
                     AA=2.0D00
                     if (I.EQ.J) AA=1.0D00
                     do 170 K=JA,JC
                        KC=IFACT(K)
                        if (I.GE.K) then
                           IK=KA+K
                        else
                           IK=0
                        end if
                        if (J.GE.K) then
                           JK=KB+K
                        else
                           JK=0
                        end if
                        do 170 L=JA,K
                           if (I.GE.L) then
                              IL=KA+L
                           else
                              IL=0
                           end if
                           if (J.GE.L) then
                              JL=KB+L
                           else
                              JL=0
                           end if
                           KL=KC+L
                           BB=2.0D00
                           if (K.EQ.L) BB=1.0D00
                           KK=KK+1
                           AJ=WJ(KK)
                           AK=WK(KK)
C
C     A  IS THE REPULSION INTEGRAL (I,J/K,L) WHERE ORBITALS I AND J ARE
C     ON ATOM II, AND ORBITALS K AND L ARE ON ATOM JJ.
C     AA AND BB ARE CORRECTION FACTORS SINCE
C     (I,J/K,L)=(J,I/K,L)=(I,J/L,K)=(J,I/L,K)
C     IJ IS THE LOCATION OF THE MATRIX ELEMENTS BETWEEN ATOMIC ORBITALS
C     I AND J.  SIMILARLY FOR IK ETC.
C
C THIS FORMS THE TWO-ELECTRON TWO-CENTER REPULSION PART OF THE FOCK
C MATRIX.  THE CODE HERE IS HARD TO FOLLOW, AND IMPOSSIBLE TO MODIFY!,
C BUT IT WORKS,
                           if (KL.LE.IJ) then
                              if (I.EQ.K.AND.AA+BB.LT.2.1D0) then
                                 BB=BB*0.5D0
                                 AA=AA*0.5D0
                                 F(IJ)=F(IJ)+BB*AJ*PTOT(KL)
                                 F(KL)=F(KL)+AA*AJ*PTOT(IJ)
                              else
                                 F(IJ)=F(IJ)+BB*AJ*PTOT(KL)
                                 F(KL)=F(KL)+AA*AJ*PTOT(IJ)
                                 A=AK*AA*BB*0.25D0
                                 F(IK)=F(IK)-A*P(JL)
                                 F(IL)=F(IL)-A*P(JK)
                                 F(JK)=F(JK)-A*P(IL)
                                 F(JL)=F(JL)-A*P(IK)
                              end if
                           end if
  170          continue
            end if
  180    continue
  190 continue
      return
C
C                    START OF MINDO/3 OPTION
C
  200 KR=0
      do 230 II=1,NUMAT
         IA=NFIRST(II)
         IB=NLAST(II)
         IM1=II-IONE
         do 220 JJ=1,IM1
            KR=KR+1
            if (LID) then
               ELREP=W(KR)
               ELEXC=ELREP
            else
               ELREP=WJ(KR)
               ELEXC=WK(KR)
            end if
            JA=NFIRST(JJ)
            JB=NLAST(JJ)
            do 210 I=IA,IB
               KA=IFACT(I)
               KK=KA+I
               do 210 K=JA,JB
                  LL=I1FACT(K)
                  IK=KA+K
                  F(KK)=F(KK)+PTOT(LL)*ELREP
                  F(LL)=F(LL)+PTOT(KK)*ELREP
  210       F(IK)=F(IK)-P(IK)*ELEXC
  220    continue
  230 continue
      return
      end
