      SUBROUTINE DCART (COORD,DXYZ)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension COORD(3,*), DXYZ(3,*)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
C***********************************************************************
C
C    DCART CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
C          CARTESIAN COORDINATES. THIS IS DONE BY FINITE DIFFERENCES.
C
C    THE MAIN ARRAYS IN DCART ARE:
C        DXYZ   ON EXIT CONTAINS THE CARTESIAN DERIVATIVES.
C
C***********************************************************************
      common /KEYWRD/ KEYWRD
      common /EULER / TVEC(3,3), ID
      common /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
      common /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
      common /DCARTC/ K1L,K2L,K3L,K1U,K2U,K3U
      common /NUMCAL/ NUMCAL
      character*241 KEYWRD
      dimension PDI(171),PADI(171),PBDI(171),
     1CDI(3,2),NDI(2),LSTOR1(6), LSTOR2(6), ENG(3)
      logical DEBUG, FORCE, MAKEP, ANADER, LARGE
      equivalence (LSTOR1(1),L1L), (LSTOR2(1), K1L)
      save CHNGE, CHNGE2, ANADER, DEBUG, FORCE
      data ICALCN/0/
      data CHNGE /1.D-4/
c
      CHNGE2=CHNGE*0.5D0
*
* CHNGE IS A MACHINE-PRECISION DEPENDENT CONSTANT
* CHNGE2=CHNGE/2
*
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
         LARGE = (INDEX(KEYWRD,'LARGE') .NE. 0)
         ANADER= (INDEX(KEYWRD,'ANALYT') .NE. 0)
         DEBUG = (INDEX(KEYWRD,'DCART') .NE. 0)
         FORCE = (INDEX(KEYWRD,'PREC')+INDEX(KEYWRD,'FORCE') .NE. 0)
      end if
      NCELLS=(L1U-L1L+1)*(L2U-L2L+1)*(L3U-L3L+1)
      do 10 I=1,6
         LSTOR2(I)=LSTOR1(I)
   10 LSTOR1(I)=0
      IOFSET=(NCELLS+1)/2
      NUMTOT=NUMAT*NCELLS
      do 20 I=1,NUMTOT
         do 20 J=1,3
   20 DXYZ(J,I)=0.D0
      if (ANADER) rewind MIOK
      KREP=0
      do 130 II=1,NUMAT
         III=NCELLS*(II-1)+IOFSET
         IM1=II
         IF=NFIRST(II)
         IM=NMIDLE(II)
         IL=NLAST(II)
         NDI(2)=NAT(II)
         do 30 I=1,3
   30    CDI(I,2)=COORD(I,II)
         do 130 JJ=1,IM1
            JJJ=NCELLS*(JJ-1)
C  FORM DIATOMIC MATRICES
            JF=NFIRST(JJ)
            JM=NMIDLE(JJ)
            JL=NLAST(JJ)
C   GET FIRST ATOM
            NDI(1)=NAT(JJ)
            MAKEP=.TRUE.
            do 120 IK=K1L,K1U
               do 120 JK=K2L,K2U
                  do 120 KL=K3L,K3U
                     JJJ=JJJ+1
                     KKK=KKK-1
                     do 40 L=1,3
   40                CDI(L,1)=COORD(L,JJ)+TVEC(L,1)*IK+TVEC(L,2)*JK
     1                                                +TVEC(L,3)*KL
                     if (.NOT. MAKEP) go to 90
                     MAKEP=.FALSE.
                     IJ=0
                     do 50 I=JF,JL
                        K=I*(I-1)/2+JF-1
                        do 50 J=JF,I
                           IJ=IJ+1
                           K=K+1
                           PADI(IJ)=PA(K)
                           PBDI(IJ)=PB(K)
   50                PDI(IJ)=P(K)
C GET SECOND ATOM FIRST ATOM INTERSECTION
                     do 80 I=IF,IL
                        L=I*(I-1)/2
                        K=L+JF-1
                        do 60 J=JF,JL
                           IJ=IJ+1
                           K=K+1
                           PADI(IJ)=PA(K)
                           PBDI(IJ)=PB(K)
   60                   PDI(IJ)=P(K)
                        K=L+IF-1
                        do 70 L=IF,I
                           K=K+1
                           IJ=IJ+1
                           PADI(IJ)=PA(K)
                           PBDI(IJ)=PB(K)
   70                   PDI(IJ)=P(K)
   80                continue
   90                continue
                     if (II.EQ.JJ) go to  120
                     if (ANADER) then
                        call ANALYT(PDI,PADI,PBDI,CDI,NDI,JF,JL,IF,IL,
     1                                                             ENG)
                        do 100 K=1,3
                           DXYZ(K,III)=DXYZ(K,III)-ENG(K)
  100                   DXYZ(K,JJJ)=DXYZ(K,JJJ)+ENG(K)
                     else
                        if ( .NOT. FORCE) then
                           CDI(1,1)=CDI(1,1)+CHNGE2
                           CDI(2,1)=CDI(2,1)+CHNGE2
                           CDI(3,1)=CDI(3,1)+CHNGE2
                           call DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM
     1                                 ,IL,NORBS,AA,1)
                        end if
                        do 110 K=1,3
                           if ( FORCE ) then
                              CDI(K,2)=CDI(K,2)-CHNGE2
                              call DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF
     1                                    ,IM,IL,NORBS,AA,1)
                           end if
                           CDI(K,2)=CDI(K,2)+CHNGE
                           call DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM
     1                                 ,IL,NORBS,EE,2)
                           CDI(K,2)=CDI(K,2)-CHNGE2
                           if ( .NOT. FORCE) CDI(K,2)=CDI(K,2)-CHNGE2
                           DERIV=(AA-EE)*23.061D0/CHNGE
                           DXYZ(K,III)=DXYZ(K,III)-DERIV
                           DXYZ(K,JJJ)=DXYZ(K,JJJ)+DERIV
  110                   continue
                     end if
  120       continue
  130 continue
      if (NNHCO.NE.0) then
C
C   NOW ADD IN MOLECULAR-MECHANICS CORRECTION TO THE H-N-C=O TORSION
C
         DEL=1.D-8
         do 160 I=1,NNHCO
            do 150 J=1,4
               do 140 K=1,3
                  COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL
                  call DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),
     1                       NHCO(4,I),ANGLE)
                  REFH=HTYPE(ITYPE)*SIN(ANGLE)**2
                  COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))+DEL*2.D0
                  call DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),
     1                       NHCO(4,I),ANGLE)
                  COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL
                  HEAT=HTYPE(ITYPE)*SIN(ANGLE)**2
                  SUM=(REFH-HEAT)/(2.D0*DEL)
                  DXYZ(K,NHCO(J,I))=DXYZ(K,NHCO(J,I))-SUM
  140          continue
  150       continue
  160    continue
      end if
      do 170 I=1,6
  170 LSTOR1(I)=LSTOR2(I)
      if (  .NOT. DEBUG) return
      write(MIOB,'(//10X,''CARTESIAN COORDINATE DERIVATIVES'',//3X,
     1 ''NUMBER  ATOM '',5X,''X'',12X,''Y'',12X,''Z'',/)')
      if (NCELLS.EQ.1) then
         write(MIOB,'(2I6,F13.6,2F13.6)')
     1         (I,NAT(I),(DXYZ(J,I),J=1,3),I=1,NUMTOT)
      else if (LARGE) then
         write(MIOB,'(2I6,F13.6,2F13.6)')
     1         (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I),J=1,3),I=1,NUMTOT)
      else
         write(MIOB,'(2I6,F13.6,2F13.6)')
     1    (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I)+DXYZ(J,I+1)+DXYZ(J,I+2)
     2                                         ,J=1,3),I=1,NUMTOT,3)
      end if
      if (ANADER) rewind MIOK
      return
      end
