      SUBROUTINE DFOCK2(F, PTOT, P, W, NUMAT, NFIRST,
     1NMIDLE, NLAST, NATI)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension F(*), PTOT(*), WJ(1), WK(1), NFIRST(*), NMIDLE(*),
     1          NLAST(*), P(*), W(*)
      double precision WJ,WK
C***********************************************************************
C
C     DFOCK2 ADDS THE 2-ELECTRON 2-CENTER REPULSION CONTRIBUTION TO
C     THE FOCK MATRIX DERIVATIVE WITHIN THE NDdo OR MINDO FORMALISMS.
C  INPUT
C     F    : 1-ELECTRON CONTRIBUTIONS DERIVATIVES.
C     PTOT : TOTAL DENSITY MATRIX.
C     P    : ALPHA OR BETA DENSITY MATRIX. = 0.5 * PTOT
C     W    : NON VANISHING TWO-ELECTRON INTEGRAL DERIVATIVES
C            (ORDERED AS DEFINED IN DHCORE).
C     NATI : # OF THE ATOM SUPPORTING THE VARYING CARTESIAN COORDINATE.
C  OUTPUT
C     F    : FOCK MATRIX DERIVATIVE WITH RESPECT TO THE CART. COORD.
C
C***********************************************************************
      common /EULER / TVEC(3,3), ID
      common /NUMCAL/ NUMCAL
      common /WORK4 / PTOT2
      common /KEYWRD/ KEYWRD
      dimension IFACT(MAXORB),
     1 I1FACT(MAXORB), JINDEX(256), KINDEX(256), IJPERM(10), LLPERM(10),
     2 PK(16), PJA(16), PJB(16), MMPERM(10),
     3 PTOT2(MAXORB,MAXORB), JJNDEX(256)
      save IFACT,I1FACT, IONE, ITYPE, LID
      character*241 KEYWRD
      logical LID
      data ITYPE /1/
      data ICALCN/0/
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
         ITYPE=0
      end if
   10 continue
      go to (20,270,70), ITYPE
   20 continue
C
C   SET UP ARRAY OF LOWER HALF TRIANGLE INDICES (PASCAL'S TRIANGLE)
C
      do 30 I=1,MAXORB
         IFACT(I)=(I*(I-1))/2
   30 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 40 I=1,4
         do 40 J=1,4
            IJ=MIN(I,J)
            JI=I+J-IJ
            do 40 K=1,4
               IK=MIN(I,K)
               KI=I+K-IK
               do 40 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
   40 JINDEX(M)=(IFACT(JI) + IJ)*10 + IFACT(LK) + KL - 10
      L=0
      do 50 I=1,4
         I1=(I-1)*4
         do 50 J=1,I
            I1=I1+1
            L=L+1
            IJPERM(L)=I1
            MMPERM(L)=IJPERM(L)-16
            LLPERM(L)=(I1-1)*16
   50 continue
      L=0
      do 60 I=1,10
         M=MMPERM(I)
         L=LLPERM(I)
         do 60 K=1,16
            L=L+1
            M=M+16
   60 JJNDEX(L)=JINDEX(M)
      LID=(ID.EQ.0)
      IONE=1
      if (ID.NE.0) IONE=0
      if (INDEX(KEYWRD,'MINDO') .NE. 0) then
         ITYPE=2
      else
         ITYPE=3
      end if
      go to 10
   70 KK=0
      L=0
      do 90 I=1,NUMAT
         IA=NFIRST(I)
         IB=NLAST(I)
         M=0
         do 80 J=IA,IB
            do 80 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)
   80    continue
   90 continue
      II=NATI
      IA=NFIRST(II)
      IB=NLAST(II)
      do 260 JJ=1,NUMAT
         if (II.EQ.JJ) go to 260
         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 100 I=1,16
                  PJA(I)=PTOT2(II,I)
  100          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
               if (IA.GT.JA) then
                  L=0
                  do 110 I=IA,IB
                     do 110 J=JA,JB
                        L=L+1
  110             PK(L)=P(IFACT(I)+J)
               else
                  L=0
                  do 120 I=IA,IB
                     do 120 J=JA,JB
                        L=L+1
  120             PK(L)=P(IFACT(J)+I)
               end if
               I1=IA
               J1=JA
               call KAB(IA,JA, PK, W(KK+1), KINDEX, F)
               IA=I1
               JA=J1
               KK=KK+100
            else if (IB-IA.GE.3) 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 140 I=0,3
                  J1=IFACT(IA+I)+IA-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
               if (IA.GT.JA) then
                  K=0
                  do 160 I=IA,IB
                     I1=IFACT(I)+JA
                     SUM=0.D0
                     do 150 J=IA,IB
                        K=K+1
                        J1=IFACT(J)+JA
  150                SUM=SUM+P(J1)*W(KK+JINDEX(K))
  160             F(I1)=F(I1)-SUM
               else
                  K=0
                  do 180 I=IA,IB
                     I1=IFACT(JA)+I
                     SUM=0.D0
                     do 170 J=IA,IB
                        K=K+1
                        J1=IFACT(JA)+J
  170                SUM=SUM+P(J1)*W(KK+JINDEX(K))
  180             F(I1)=F(I1)-SUM
               end if
               KK=KK+10
            else if (JB-JA.GE.3) 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 200 I=0,3
                  J1=IFACT(JA+I)+JA-1
                  do 190 J=0,I-1
                     K=K+1
                     J1=J1+1
                     F(J1)=F(J1)+PTOT(LL)*W(KK+K)
  190             SUMOFF=SUMOFF+PTOT(J1)*W(KK+K)
                  J1=J1+1
                  K=K+1
                  F(J1)=F(J1)+PTOT(LL)*W(KK+K)
  200          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
               if (IA.GT.JA) then
                  K=IFACT(IA)+JA
                  J=0
                  do 220 I=K,K+3
                     SUM=0.D0
                     do 210 L=K,K+3
                        J=J+1
  210                SUM=SUM+P(L)*W(KK+JINDEX(J))
  220             F(I)=F(I)-SUM
               else
                  J=0
                  do 240 K=JA,JA+3
                     I=IFACT(K)+IA
                     SUM=0.D0
                     do 230 LL=JA,JA+3
                        L=IFACT(LL)+IA
                        J=J+1
  230                SUM=SUM+P(L)*W(KK+JINDEX(J))
  240             F(I)=F(I)-SUM
               end if
               KK=KK+10
            else
C
C                         LIGHT-ATOM - LIGHT-ATOM
C
               I1=I1FACT(IA)
               J1=I1FACT(JA)
               F(I1)=F(I1)+PTOT(J1)*W(KK+1)
               F(J1)=F(J1)+PTOT(I1)*W(KK+1)
               if (IA.GT.JA) then
                  IJ=I1+JA-IA
                  F(IJ)=F(IJ)-P   (IJ)*W(KK+1)
               else
                  IJ=J1+IA-JA
                  F(IJ)=F(IJ)-P   (IJ)*W(KK+1)
               end if
               KK=KK+1
            end if
         else
            DREP=WJ(KK+1)
            do 250 I=IA,IB
               KA=IFACT(I)
               do 250 J=IA,I
                  KB=IFACT(J)
                  IJ=KA+J
                  AA=2.0D00
                  if (I.EQ.J) AA=1.0D00
                  do 250 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 250 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
  250       continue
         end if
  260 continue
C
      return
  270 KR=0
      II=NATI
      IA=NFIRST(II)
      IB=NLAST(II)
      do 290 JJ=1,NUMAT
         if (JJ.EQ.II) go to 290
         KR=KR+1
         ELREP=W(KR)
         JA=NFIRST(JJ)
         JB=NLAST(JJ)
         do 280 I=IA,IB
            KA=IFACT(I)
            KK=KA+I
            do 280 K=JA,JB
               LL=I1FACT(K)
               if (JA.LT.IA) then
                  IK=KA+K
               else
                  IK=LL+I-K
               end if
               F(KK)=F(KK)+PTOT(LL)*ELREP
               F(LL)=F(LL)+PTOT(KK)*ELREP
  280    F(IK)=F(IK)-P(IK)*ELREP
  290 continue
      return
      end
