C
C $Id: denful.F,v 1.3 1998/07/16 16:39:33 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE DENFUL(NORBS,EVEC1)
C
C     COMPUTES THE GLOBAL DENSITY MATRIX (PDIAG,PDIAT) FOR THE FULL SYSTEM
C     IN A CALCULATION THAT DOES NOT INVOLVE SEPARATE SUBSYSTEMS.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      DIMENSION EVEC1(NORBS,*)
C
      NOCC = NELEC/2
      NHALF = NORBS/2
      IF(NOCC.GT.NHALF)THEN
C
C       DENSITY MATRIX WILL BE FORMED BY SUBTRACTING ELECTRONS FROM
C       A SATURATED SYSTEM.  THIS IS USUALLY THE CASE.
C
        NOCC1 = NOCC + 1
        NOCC2 = NORBS
        PSIGN = -2.0D0
      ELSE
C
C       DENSITY MATRIX WILL BE FORMED BY ADDING ELECTRONS TO AN EMPTY SYSTEM.
C
        NOCC1 = 1
        NOCC2 = NOCC
        PSIGN = 2.0
      ENDIF
      IIMAX = IIMAT(NATOMS+1)-1
      IPMAX = IP1(NATOMS+1)
      IJMAX = IJMAT(IPMAX)-1
      DO 10 I=1,IIMAX
        PDIAG(I) = 0.0D0
 10   CONTINUE
      DO 20 IJ=1,IJMAX
        PDIAT(IJ) = 0.0D0
 20   CONTINUE
      IF(PSIGN.LT.0.0D0)THEN
C
C       ASSIGN A 'SATURATED' DENSITY MATRIX WITH EACH DIAGONAL=2.0.
C
        II = 0
        DO 50 I=1,NATOMS
          NORBSI = NATORB(IATNUM(I))
C
C         SKIP IF DUMMY OR SPARKLE.
C
          IF(NORBSI.EQ.0) GO TO 50
          DO 40 IDIAG=1,NORBSI
            II = II + IDIAG
            PDIAG(II) = 2.0D0
 40       CONTINUE
 50     CONTINUE
      ENDIF
      NATMS = IATOM1(2)-1
C
C     DIAGONAL BLOCKS FIRST.
C
      DO 100 I=1,NATMS
        IATM = IATOMS(I)
        NORBSI = NATORB(IATNUM(IATM))
        IJP = IIMAT(IATM)
        I0 = IORB1(I)-1
        J0 = I0
        DO 80 IP=1,NORBSI
          IEVEC1 = I0+IP
          DO 70 JP=1,IP
            PIJ = 0.0D0
            JEVEC1 = J0+JP
            DO 60 L=NOCC1,NOCC2
              PIJ = PIJ + EVEC1(IEVEC1,L)*EVEC1(JEVEC1,L)
 60         CONTINUE
            PDIAG(IJP) = PDIAG(IJP) + PSIGN*PIJ
            IJP = IJP + 1
 70       CONTINUE
 80     CONTINUE
 100  CONTINUE
C
C     NOW OFF-DIAGONAL BLOCKS.
C
      IF(NATMS.EQ.1) RETURN
      NPAIRS = IP1(NATOMS+1)-1
      IATM1 = IATOMS(IATOM1(1))
      DO 200 I=2,NATMS
        IATM = IATOMS(I)
        NORBSI = NATORB(IATNUM(IATM))
        I0 = IORB1(I)-1
        DO 180 J=1,I-1
          JATM = IATOMS(J)
          NORBSJ = NATORB(IATNUM(JATM))
C
C         FIND POSTION IN PAIRLIST OF (IATM,JATM) PAIR.
C
          CALL IJFIND(NPAIRS,IATM,JATM,IJADDR)
          IF(IJADDR.EQ.0) GO TO 180
          IJPMAT = IJMAT(IJADDR)
          J0 = IORB1(J)-1
          DO 160 IP=1,NORBSI
            IEVEC1 = I0+IP
            DO 150 JP=1,NORBSJ
              PIJ = 0.0D0
              JEVEC1 = J0+JP
              DO 140 L=NOCC1,NOCC2
                PIJ = PIJ + EVEC1(IEVEC1,L)*EVEC1(JEVEC1,L)
 140          CONTINUE
              PDIAT(IJPMAT) = PDIAT(IJPMAT) + PSIGN*PIJ
              IJPMAT = IJPMAT + 1
 150        CONTINUE
 160      CONTINUE
 180    CONTINUE
 200  CONTINUE
      RETURN
      END