      SUBROUTINE HCORE (COORD,H,W, WJ,WK,ENUCLR)
      implicit double precision (A-H,O-Z)
      logical FLDON
      include 'SIZES'
      include 'MOP6FILE'
      dimension COORD(3,*),H(*), WJ(N2ELEC), WK(N2ELEC), W(N2ELEC)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
     3       /MOLORB/ USPD(MAXORB),DUMY(MAXORB)
     4       /KEYWRD/ KEYWRD
      common /EULER / TVEC(3,3), ID
      common /MULTIP/ DD(107),QQ(107),AM(107),AD(107),AQ(107)
      common /CORE  / CORE(107)
      common /FIELD / EFIELD(3)
      common /NUMCAL/ NUMCAL
************************************************************************
C
C   HCORE GENERATES THE ONE-ELECTRON MATRIX AND TWO ELECTRON INTEGRALS
C         FOR A GIVEN MOLECULE WHOSE GEOMETRY IS DEFINED IN CARTESIAN
C         COORDINATES.
C
C  ON INPUT  COORD   = COORDINATES OF THE MOLECULE.
C
C  ON OUTPUT  H      = ONE-ELECTRON MATRIX.
C             W      = TWO-ELECTRON INTEGRALS.
C             ENUCLR = NUCLEAR ENERGY
************************************************************************
      character*241 KEYWRD
      logical FIRST,DEBUG
      save FIRST, IONE, CUTOFF, DEBUG
      dimension E1B(10),E2A(10),DI(9,9), WJD(100), WKD(100)
      data ICALCN/0/
      FIRST=(ICALCN.NE.NUMCAL)
      ICALCN=NUMCAL
      if (FIRST) then
         IONE=1
         CUTOFF=1.D10
         if (ID.NE.0) CUTOFF=60.D0
         if (ID.NE.0) IONE=0
         DEBUG=(INDEX(KEYWRD,'HCORE') .NE. 0)
         EFIELD(1)=0.0d0  !
         EFIELD(2)=0.0d0  !  according to errata of May 1991
         EFIELD(3)=0.0d0  !
      end if
      FLDON = .FALSE.
      if ((EFIELD(1).NE.0.0D00) .OR.
     1    (EFIELD(2).NE.0.0D00) .OR.
     2    (EFIELD(3).NE.0.0D00))       then
         FLDCON = 51.4257D00
         FLDON = .TRUE.
      end if
      do 10 I=1,(NORBS*(NORBS+1))/2
   10 H(I)=0.D0
      ENUCLR=0.D0
      KR=1
      do 110 I=1,NUMAT
         IA=NFIRST(I)
         IB=NLAST(I)
         IC=NMIDLE(I)
         NI=NAT(I)
C
C FIRST WE FILL THE DIAGONALS, AND OFF-DIAGONALS ON THE SAME ATOM
C
         do 30 I1=IA,IB
            I2=I1*(I1-1)/2+IA-1
            do 20 J1=IA,I1
               I2=I2+1
               H(I2)=0.D0
               if (FLDON) then
                  IO1 = I1 - IA
                  JO1 = J1 - IA
                  if ((JO1.EQ.0).AND.(IO1.EQ.1)) then
                     HTERME = -0.529177D00*DD(NI)*EFIELD(1)*FLDCON
                     H(I2) = HTERME
                  end if
                  if ((JO1.EQ.0).AND.(IO1.EQ.2)) then
                     HTERME = -0.529177D00*DD(NI)*EFIELD(2)*FLDCON
                     H(I2) = HTERME
                  end if
                  if ((JO1.EQ.0).AND.(IO1.EQ.3)) then
                     HTERME = -0.529177D00*DD(NI)*EFIELD(3)*FLDCON
                     H(I2) = HTERME
                  end if
               end if
   20       continue
            H(I2) = USPD(I1)
            if (FLDON) then
               FNUC = -(EFIELD(1)*COORD(1,I) + EFIELD(2)*COORD(2,I) +
     1                  EFIELD(3)*COORD(3,I))*FLDCON
               H(I2) = H(I2) + FNUC
            end if
   30    continue
C
C   FILL THE ATOM-OTHER ATOM ONE-ELECTRON MATRIX<PSI(LAMBDA)|PSI(SIGMA)>
C
         IM1=I-IONE
         do 100 J=1,IM1
            HALF=1.D0
            if (I.EQ.J) HALF=0.5D0
            JA=NFIRST(J)
            JB=NLAST(J)
            JC=NMIDLE(J)
            NJ=NAT(J)
            call H1ELEC(NI,NJ,COORD(1,I),COORD(1,J),DI)
            I2=0
            do 40 I1=IA,IB
               II=I1*(I1-1)/2+JA-1
               I2=I2+1
               J2=0
               JJ=MIN(I1,JB)
               do 40 J1=JA,JJ
                  II=II+1
                  J2=J2+1
   40       H(II)=H(II)+DI(I2,J2)
C
C   CALCULATE THE TWO-ELECTRON INTEGRALS, W; THE ELECTRON NUCLEAR TERMS
C   E1B AND E2A; AND THE NUCLEAR-NUCLEAR TERM ENUC.
C
            if (ID.EQ.0) then
               call ROTATE(NI,NJ,COORD(1,I),COORD(1,J),
     1                     W(KR), KR,E1B,E2A,ENUC,CUTOFF)
            else
               KRO=KR
               call SOLROT(NI,NJ,COORD(1,I),COORD(1,J),
     1                     WJD, WKD,KR,E1B,E2A,ENUC,CUTOFF)
               JJ=0
               do 50 II=KRO,KR-1
                  JJ=JJ+1
                  WJ(II)=WJD(JJ)
   50          WK(II)=WKD(JJ)
            end if
            ENUCLR = ENUCLR + ENUC
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM I.
C
            I2=0
            do 60 I1=IA,IC
               II=I1*(I1-1)/2+IA-1
               do 60 J1=IA,I1
                  II=II+1
                  I2=I2+1
   60       H(II)=H(II)+E1B(I2)*HALF
            do  70 I1=IC+1,IB
               II=(I1*(I1+1))/2
   70       H(II)=H(II)+E1B(1)*HALF
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM J.
C
            I2=0
            do 80 I1=JA,JC
               II=I1*(I1-1)/2+JA-1
               do 80 J1=JA,I1
                  II=II+1
                  I2=I2+1
   80       H(II)=H(II)+E2A(I2)*HALF
            do 90 I1=JC+1,JB
               II=(I1*(I1+1))/2
   90       H(II)=H(II)+E2A(1)*HALF
  100    continue
  110 continue
      if ( .NOT. DEBUG) return
      write(MIOB,'(//10X,''ONE-ELECTRON MATRIX FROM HCORE'')')
      call VECPRT(H,NORBS)
      J=MIN(400,KR)
      if (ID.EQ.0) then
         write(MIOB,'(//10X,''TWO-ELECTRON MATRIX IN HCORE''/)')
         write(MIOB,120) (W(I),I=1,J)
      else
         write(MIOB,'(//10X,''TWO-ELECTRON J MATRIX IN HCORE''/)')
         write(MIOB,120) (WJ(I),I=1,J)
         write(MIOB,'(//10X,''TWO-ELECTRON K MATRIX IN HCORE''/)')
         write(MIOB,120) (WK(I),I=1,J)
      end if
  120 format(10F8.4)
      return
      end
