C
C $Id: densub.F,v 1.4 1998/07/16 16:39:33 jjv5 Exp $
C
C------------------------------------------------------------------------
      SUBROUTINE DENSUBfdm
C
C     ADDS IN DENSITY MATRIX CONTRIBUTION FROM SUBSYSTEM K.  IT IS
C     ASSUMED THAT IT WILL BE FASTER TO SUBTRACT ELECTRONS FROM A
C     SATURATED DENSITY MATRIX, WHICH IS JUST 2.0 TIMES THE IDENTITY
C     MATRIX.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
C     LOCAL:
C
      DIMENSION HOLE(MSORB)
C
C     INITIALIZE A SATURATED DENSITY MATRIX IF THIS IS THE FIRST
C     SUBSYSTEM.
C

      IPMAX = IP1(NATOMS+1)
      IJMAX = IJMAT(IPMAX)-1
      DO 30 I=1,NATOMS
         NORBSI = NATORB(IATNUM(I))
         IF(NORBSI.EQ.0) GO TO 30
         IJ = IIMAT(I)
         DO 20 II=1,NORBSI
            DO 10 JJ=1,II
               PDIAG(IJ) = 0.0D0
               IF((JJ.EQ.II) .AND. (myid .EQ. 0)) PDIAG(IJ) = 2.0D0
               IJ = IJ+1
 10         CONTINUE
 20      CONTINUE
 30   CONTINUE
      DO 40 IJ=1,IJMAX
         PDIAT(IJ) = 0.0D0
 40   CONTINUE

C
C     DETERMINE LOWER MO LIMIT FOR ELECTRON SUBTRACTION AND ASSIGN
C     HOLE OCCUPATION NUMBERS FROM FERMI OCCUPATION NUMBERS.
C
      k0 = 0
      do k=1,nsub
         NORBSK = IORBPT(K+1)-IORBPT(K)
         NLOWER = 0
         DO 50 I=IORBPT(K),IORBPT(K+1)-1
            NLOWER = NLOWER + 1
            IF(FERMI(I).LT.0.999999D0) GO TO 100
 50      CONTINUE
 100     L0 = IORBPT(K)-1
         DO 110 L=NLOWER,NORBSK
            HOLE(L) = FERMI(L0+L) - 1.0D0
 110     CONTINUE
C
C     DIAGONAL BLOCKS FIRST.
C
         I1 = IATOM1(K)
         I2 = IATOM1(K+1)-1
         DO 200 I=I1,I2
C
C       SKIP DENSITY CONTRIBUTION IF THIS IS A BUFFER ATOM.
C
            IF(IABUFF(I).NE.0) GO TO 200
C
            IATM = IATOMS(I)
            NORBSI = NATORB(IATNUM(IATM))
            NSUBSI = ISUB1(IATM+1) - ISUB1(IATM)
            DII = 2.0D0/NSUBSI
            IJP = IIMAT(IATM)
            I0 = IORB1(I)-1
            J0 = I0
            DO 170 IP=1,NORBSI
               IEVEC = I0+IP
               ig0 = ievec + k0
               DO 160 JP=1,IP
                  JEVEC = J0+JP
                  jg0 = jevec + k0
                  PIJ = 0.0D0
                  DO 150 L=NLOWER,NORBSK
                     il = (l-1)*norbsk
                     PIJ = PIJ + EVEC(il+ig0)*EVEC(il+jg0)*HOLE(L)
 150              CONTINUE
                  PDIAG(IJP) = PDIAG(IJP) + PIJ*DII
                  IJP = IJP + 1
 160           CONTINUE
 170        CONTINUE
 200     CONTINUE
C
C     NOW OFF-DIAGONAL BLOCKS.
C
         NATMS = IATOM1(K+1) - IATOM1(K)
         IF(NATMS.EQ.1) RETURN
         NPAIRS = IP1(NATOMS+1)-1
         IATM0 = IATOM1(K)-1
         DO 300 I=2,NATMS
            I0I = IATM0+I
C
C     NO BONDING WITH I0I IF IT'S AN OUTER LAYER BUFFER ATOM.
C
            IF(IABUFF(I0I).EQ.2) GO TO 300
C
            IATM = IATOMS(I0I)
            NORBSI = NATORB(IATNUM(IATM))
            I0 = IORB1(I0I)-1
            DO 280 J=1,I-1
               I0J = IATM0+J
C
C     NO BONDING WITH I0J IF IT'S AN OUTER LAYER BUFFER ATOM OR
C     IF BOTH I0J AND I0I HAVE ANY SORT OF BUFFER STATUS.
C
               IJSUM = IABUFF(I0I) + IABUFF(I0J)
               IF(IJSUM.GE.2) GO TO 280
C
               JATM = IATOMS(I0J)
C
C         FIND POSTION IN PAIRLIST OF (IATM,JATM) PAIR.
C
               CALL IJFIND(NPAIRS,IATM,JATM,IJADDR)
C
C     SKIP THIS PAIR IF IT HASN'T BEEN STORED AS A BONDED PAIR.
C
               IF(IJADDR.EQ.0) GO TO 280
C
               NORBSJ = NATORB(IATNUM(JATM))
               J0 = IORB1(I0J)-1
               IJPMAT = IJMAT(IJADDR)
               DIJ = 2.0D0/NSHARE(1,IJADDR)
               DO 260 IP=1,NORBSI
                  IEVEC = I0+IP
                  ig0 = ievec + k0
                  DO 250 JP=1,NORBSJ
                     PIJ = 0.0D0
                     JEVEC = J0+JP
                     jg0 = jevec + k0
                     DO 240 L=NLOWER,NORBSK
                        il = (l-1)*norbsk
                        PIJ = PIJ + EVEC(il+ig0)*EVEC(il+jg0)*HOLE(L)
 240                 CONTINUE
                     PDIAT(IJPMAT) = PDIAT(IJPMAT) + PIJ*DIJ
                     IJPMAT = IJPMAT + 1
 250              CONTINUE
 260           CONTINUE
 280        CONTINUE
 300     CONTINUE
         k0 = k0 + norbsk*norbsk
      enddo
      RETURN
      END