*Deck EFill
      SUBROUTINE EFILL(ISHELL,JSHELL,KSHELL,LSHELL,LA,LB,LC,LD,AX1,
     $    ISCF,DM,DN,E,DMAX,NBasis,NonCor,NonC1,FHam,AJ,AK,CMO,INDDO,
     $    AtmChg,FirstS)
      Implicit Real*8(A-H,O-Z)
C
C     ROUTINE TO PLUCK DENSITY MATRIX CONTRIBUTIONS ACCORDING
C     TO SHELL NUMBERS (ISHELL,...,LSHELL) AND STORE INTO LOCAL
C     ARRAY E.  ALSO, THE MAX DENSITY MATRIX ELEMENT FOR THIS BLOCK
C     IS DETERMINED.
C
C     ARGUMENTS:
C
C     ISHELL TO LSHELL ... SHELL NUMBERS OF THE FOUR SHELLS.
C     ISCF             ... SCF MODE FLAG.  SEE COMMENTS IN
C                          D2ESP FOR DETAILS.
C     DM AND DN        ... ARRAYS CONTAINING DENSITY MATRICES.
C                          EXACTLY WHAT THEY CONTAIN DEPENDS ON
C                          ISCF.  REFER TO D2ESP FOR DETAILS.
C     E                ... ARRAY OF DIMENSION 256, FILLED WITH
C                          COMBINATIONS OF DENSITY MATRIX ELEMENTS.
C     DMAX             ... MAX COMBINATION FORMED.  USEFUL LATER ON
C                          WHEN DECIDING JUST WHAT TO KEEP.
C     NBasis,etc       ... only used for GVB wavefunctions ... the J
C                          and K coefficients are transformed to the
C                          AO basis.
C     INDDO,AtmChg,FirstS used for NDDO gradients.
C
      INTEGER SHELLA,SHELLN,SHELLT,SHELLC,AOS,AON
      LOGICAL OPEN,COMPLX,FirstS(1),FrstKL
      DIMENSION DM(1),DN(1),E(1),FHam(NonC1),AJ(NonC1,NonC1),
     $    AK(NonC1,NonC1),CMO(NBasis,NonC1),AtmChg(1)
      COMMON/B/EXX(1200),C1(1200),C2(1200),C3(1200),
     $         X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400),
     $         SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP
      DIMENSION C4(400),SHLADF(400)
      EQUIVALENCE(C4(1),C3(401)),(SHLADF(1),C3(801))
      Data Zero,P25,H/0.D0,0.25D0,0.5D0/,One/1.0d0/,Two/2.0d0/,
     $     Four/4.0d0/
C
C     SET UP LOCAL CONTROL VARIABLES FOR WITHIN LOOPS.
      ISCFP=ISCF+1
      OPEN=ISCF.EQ.1
      COMPLX=ISCF.EQ.2
C
C     DETERMINE STARTING ADDRESSES OF BASIS FUNCTIONS WITHIN SHELLS.
      IAOS=AOS(ISHELL)-1
      JAOS=AOS(JSHELL)-1
      KAOS=AOS(KSHELL)-1
      LAOS=AOS(LSHELL)-1
C
C     LOOP OVER ALL BASIS FUNCTIONS FOR THIS SHELL SET, AND STORE
C     DENSITY MATRIX CONTRIBUTIONS INTO E.
      DMAX=ZERO
      D12B=Zero
      D34B=Zero
      DO 20 L=1,LD
      LAS=LAOS+L
      LLL=LAS*(LAS-1)/2
C
      DO 20 K=1,LC
      KAS=KAOS+K
      LKK=KAS*(KAS-1)/2
      LKL=LKK+LAS
      IF(KAS.LT.LAS) LKL=LLL+KAS
      FrstKL=FirstS(KAS).and.FirstS(LAS)
      D34=DM(LKL)
      IF(OPEN) D34B=DN(LKL)
C
      DO 20 J=1,LB
      ID=16*J+4*K+L-84
      JAS=JAOS+J
      LJJ=JAS*(JAS-1)/2
      LJK=LJJ+KAS
      IF(JAS.LT.KAS) LJK=LKK+JAS
      D23=DM(LJK)
      LJL=LJJ+LAS
      IF(JAS.LT.LAS) LJL=LLL+JAS
      D24=DM(LJL)
      IF(.NOT.(OPEN.OR.COMPLX)) GO TO 10
      D23B=DN(LJK)
      D24B=DN(LJL)
      IF(.NOT.COMPLX) GO TO 10
      IF(KAS.LT.JAS) D23B=-D23B
      IF(JAS.LT.LAS) D24B=-D24B
   10 CONTINUE
C
      DO 20 I=1,LA
      ID=ID+64
      IAS=IAOS+I
      LII=IAS*(IAS-1)/2
      LIJ=LII+JAS
      IF(IAS.LT.JAS) LIJ=LJJ+IAS
      D12=DM(LIJ)
      LIK=LII+KAS
      IF(IAS.LT.KAS) LIK=LKK+IAS
      D13=DM(LIK)
      LIL=LII+LAS
      IF(IAS.LT.LAS) LIL=LLL+IAS
      D14=DM(LIL)
C
C     BRANCH TO PARTICULAR CODE BASED ON ISCF.
      GO TO(1,2,3,4,5),ISCFP
C
C     REAL-RHF, CLOSED SHELL.
 1    D1234=(D12*D34-P25*(D13*D24+D23*D14))*AX1
      GO TO 15
C
C     REAL-UHF.
 2    D12B=DN(LIJ)
      D13B=DN(LIK)
      D14B=DN(LIL)
      D1234=((D12+D12B)*(D34+D34B)
     *  -H*(D13*D24+D23*D14+D13B*D24B+D23B*D14B))*AX1
      GO TO 15
C
C     COMPLEX-RHF, CLOSED SHELL.
 3    D14B=DN(LIL)
      D13B=DN(LIK)
      IF(IAS.LT.LAS) D14B=-D14B
      IF(KAS.LT.IAS) D13B=-D13B
      D1234=(D12*D34-P25*(D13*D24+D23*D14-D13B*D24B-D23B*D14B))*AX1
      GO TO 15
C
C     Complex UHF ... not here yet.
 4    Goto 15
C
C     GVB.
    5 D1234 = Four*D12*D34 - D13*D24 - D14*D23
      Do 100 IOrb = 1, NonCor
  100   D1234 = D1234 + 
     $    FHam(IOrb)*(Four*(D12*CMO(KAS,IOrb)*CMO(LAS,IOrb)+
     $    D34*CMO(IAS,IOrb)*CMO(JAS,IOrb))-
     $    D13*CMO(JAS,IOrb)*CMO(LAS,IOrb)-
     $    D14*CMO(JAS,IOrb)*CMO(KAS,IOrb)-
     $    D24*CMO(IAS,IOrb)*CMO(KAS,IOrb)-
     $    D23*CMO(IAS,IOrb)*CMO(LAS,IOrb))
      Do 110 IOrb = 1, NonCor
        Do 110 JOrb = 1, NonCor
  110     D1234 = D1234+Two*CMO(IAS,IOrb)*CMO(JAS,IOrb)*CMO(KAS,JOrb)*
     $      CMO(LAS,JOrb)*AJ(IOrb,JOrb)+CMO(IAS,IOrb)*CMO(JAS,JOrb)*
     $      CMO(KAS,IOrb)*CMO(LAS,JOrb)*AK(IOrb,JOrb)+CMO(IAS,IOrb)*
     $      CMO(JAS,JOrb)*CMO(KAS,JOrb)*CMO(LAS,IOrb)*AK(IOrb,JOrb)
      D1234=AX1*D1234
C
C     STORE COMPUTED VALUE INTO E, AND CHECK DMAX.
   15 D1234 = D1234
      If(FrstKL) D1234 = D1234 - 
     $    AtmChg(JAn(KShell))*(D12+D12B)*AX1
      If(FirstS(IAS).and.FirstS(JAS)) D1234 = D1234 -
     $    AtmChg(JAn(IShell))*(D34+D34B)*AX1
      E(ID)=D1234
      D1234=DABS(D1234)
      IF(D1234.GT.DMAX) DMAX=D1234
   20 CONTINUE
C
C     ALL DONE, RETURN.
      RETURN
      END
