C
C $Id: densub.F,v 1.4 1998/07/16 16:39:33 jjv5 Exp $
C
C------------------------------------------------------------------------
      SUBROUTINE DENSUB(K,NORBSK,EVEC1)
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"
#ifdef MPI
#include "mpif.h"
#endif
      DIMENSION EVEC1(NORBSK,*)
C
C     LOCAL:
C
      DIMENSION HOLE(MSORB)
      logical dodens
C
C     INITIALIZE A SATURATED DENSITY MATRIX IF THIS IS THE FIRST
C     SUBSYSTEM.
C

      if (.not.fdmx) then
         dodens = .true.
      else
         if (notfroz(k)) then
            dodens = .true.
         else
            dodens = .false.
         endif
      endif

#ifdef MPI 
C     JV First time through each PE must init density matrix
      if (K .eq. my_subs(1) ) then
#else
      IF (K.EQ.1) THEN
#endif

            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
C     JV Only master inits diag so that we don't sum from all PE's
#ifdef MPI 
                     IF((JJ.EQ.II) .AND. (myid .EQ. 0))
#else
                     IF((JJ.EQ.II))
#endif
     $                    PDIAG(IJ) = 2.0D0
                     IJ = IJ+1
 10               CONTINUE
 20            CONTINUE
 30         CONTINUE
            DO 40 IJ=1,IJMAX
               PDIAT(IJ) = 0.0D0
 40         CONTINUE
      ENDIF

      if (dodens) then
C     .  calculate the density matrix,
C     .  if fdmx is true, this can be skipped: density matrix is held constant
C
C
C     DETERMINE LOWER MO LIMIT FOR ELECTRON SUBTRACTION AND ASSIGN
C     HOLE OCCUPATION NUMBERS FROM FERMI OCCUPATION NUMBERS.
C
         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
               IEVEC1 = I0+IP
               DO 160 JP=1,IP
                  JEVEC1 = J0+JP
                  PIJ = 0.0D0
                  DO 150 L=NLOWER,NORBSK
                     PIJ = PIJ + EVEC1(IEVEC1,L)*EVEC1(JEVEC1,L)*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
                  IEVEC1 = I0+IP
                  DO 250 JP=1,NORBSJ
                     PIJ = 0.0D0
                     JEVEC1 = J0+JP
                     DO 240 L=NLOWER,NORBSK
                        PIJ = PIJ + EVEC1(IEVEC1,L)
     $                       *EVEC1(JEVEC1,L)*HOLE(L)
 240                 CONTINUE
                     PDIAT(IJPMAT) = PDIAT(IJPMAT) + PIJ*DIJ
                     IJPMAT = IJPMAT + 1
 250              CONTINUE
 260           CONTINUE
 280        CONTINUE
 300     CONTINUE
      endif

C     just in case fdmx is true: copy the old density matrix elements
      if (fdmx.and.(k.eq.1)) then
         i=1
 400     if (i.lt.nfdiag) then
            j1 = ifdiag(i)
            j2 = ifdiag(i+1)
            do j=j1,j2
               pdiag(j) = piiold(j)
            enddo
            i = i + 2
            goto 400
         endif
         i=1
 410     if (i.lt.nfdiat) then
            j1 = ifdiat(i)
            j2 = ifdiat(i+1)
            do j=j1,j2
               pdiat(j) = pijold(j)
            enddo
            i = i + 2
            goto 410
         endif
      endif
      RETURN
      END