      SUBROUTINE DHCORE (COORD,H,W,ENUCLR,NATI,NATX,STEP)
      implicit double precision  (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension COORD(3,*),H(*),W(*)
C
C  DHCORE GENERATES THE 1-ELECTRON  AND 2-ELECTRON INTEGRALS DERIVATIVES
C         WITH RESPECT TO THE CARTESIAN COORDINATE COORD (NATX,NATI).
C
C  INPUT
C      COORD     : CARTESIAN  COORDINATES OF THE MOLECULE.
C      NATI,NATX : INDICES OF THE MOVING COORDINATE.
C      STEP      : STEP SIZE OF THE 2-POINTS FINITE DifFERENCE.
C  OUTPUT
C      H         : 1-ELECTRON INTEGRALS DERIVATIVES (PACKED CANONICAL).
C      W         : 2-ELECTRON INTEGRALS DERIVATIVES (ORDERED AS REQUIRED
C                             IN DFOCK2 AND DIJKL1).
C      ENUCLR    : NUCLEAR ENERGY DERIVATIVE.
C
      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
     5       /WMATRX/ WDUMMY(N2ELEC*2)
      character*241 KEYWRD
      logical FIRST,MINdo
      dimension E1B(10),DE1B(10),E2A(10),DE2A(10)
     1         ,DI(9,9),DDI(9,9),WJD(101),DWJD(101)
     2,NB(0:8)
      data NB/1,0,0,10,0,0,0,0,45/
      data FIRST/.TRUE./
      if (FIRST) then
         IONE=1
         CUTOFF=1.D10
         FIRST=.FALSE.
         MINDO=INDEX(KEYWRD,'MINDO') .NE. 0
      end if
      do 10 I=1,(NORBS*(NORBS+1))/2
   10 H(I)=0
      ENUCLR=0.D0
      KR=1
      NROW=0
      I=NATI
      CSAVE=COORD(NATX,NATI)
      IA=NFIRST(NATI)
      IB=NLAST(NATI)
      IC=NMIDLE(NATI)
      NI=NAT(NATI)
      NROW=-NB(IB-IA)
      do 20 J=1,NUMAT
   20 NROW=NROW+NB(NLAST(J)-NFIRST(J))
      NCOL=NB(NLAST(NATI)-NFIRST(NATI))
      NBAND2=0
      do 120 J=1,NUMAT
         if (J.EQ.NATI) go to 120
         JA=NFIRST(J)
         JB=NLAST(J)
         JC=NMIDLE(J)
         NJ=NAT(J)
         COORD(NATX,NATI)=CSAVE+STEP
         call H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DI)
C
C  THE FOLLOWING STYLE WAS NECESSARY TO GET ROUND A BUG IN THE
C  GOULD COMPILER
C
         COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
         call H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DDI)
C
C     FILL THE ATOM-OTHER ATOM ONE-ELECTRON MATRIX.
C
         I2=0
         if (IA.GT.JA) then
            do 30 I1=IA,IB
               IJ=I1*(I1-1)/2+JA-1
               I2=I2+1
               J2=0
               do 30 J1=JA,JB
                  IJ=IJ+1
                  J2=J2+1
   30       H(IJ)=H(IJ)+(DI(I2,J2)-DDI(I2,J2))
         else
            do 40 I1=JA,JB
               IJ=I1*(I1-1)/2+IA-1
               I2=I2+1
               J2=0
               do 40 J1=IA,IB
                  IJ=IJ+1
                  J2=J2+1
   40       H(IJ)=H(IJ)+(DI(J2,I2)-DDI(J2,I2))
         end if
C
C     CALCULATE THE TWO-ELECTRON INTEGRALS, W; THE ELECTRON NUCLEAR TERM
C     E1B AND E2A; AND THE NUCLEAR-NUCLEAR TERM ENUC.
C
         KRO=KR
         NBAND1=NBAND2+1
         NBAND2=NBAND2+NB(NLAST(J)-NFIRST(J))
         if (MINDO) then
            COORD(NATX,NATI)=CSAVE+STEP
            call ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
            KR=KRO
            COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
            call ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
            if (KR.GT.KRO) then
               do 50 K=1,KR-KRO+1
   50          W(KRO+K-1)=WJD(K)-DWJD(K)
            end if
         else
            COORD(NATX,NATI)=CSAVE+STEP
            call ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
            KR=KRO
            COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
            call ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
     1               ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
            if (KR.GT.KRO) then
               do 60 K=1,KR-KRO+1
   60          WJD(K)=WJD(K)-DWJD(K)
               J7=0
               do 70 I1=KRO,KR
                  J7=J7+1
   70          W(I1)=WJD(J7)
            end if
         end if
         COORD(NATX,NATI)=CSAVE
         ENUCLR = ENUCLR + ENUC-DENUC
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM I.
C
         I2=0
         do 80 I1=IA,IC
            II=I1*(I1-1)/2+IA-1
            do 80 J1=IA,I1
               II=II+1
               I2=I2+1
   80    H(II)=H(II)+E1B(I2)-DE1B(I2)
C     CONTRIB D, CNDO.
         do 90 I1=IC+1,IB
            II=(I1*(I1+1))/2
   90    H(II)=H(II)+E1B(1)-DE1B(1)
C
C   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM J.
C
         I2=0
         do 100 I1=JA,JC
            II=I1*(I1-1)/2+JA-1
            do 100 J1=JA,I1
               II=II+1
               I2=I2+1
  100    H(II)=H(II)+E2A(I2)-DE2A(I2)
C     CONTRIB D, CNDO.
         do 110 I1=JC+1,JB
            II=(I1*(I1+1))/2
  110    H(II)=H(II)+E2A(1)-DE2A(1)
  120 continue
C
C   'SIZE' OF H IS NROW * NCOL
C
      return
      end
