C
C $Id: gcartmc.F,v 1.3 1998/07/16 16:39:51 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE GCARTMC(GXYZ)
C
C     SUBROUTINE TO COMPUTE CARTESIAN ENERGY GRADIENT USING VARIATIONAL
C     FINITE-DIFFERENCE DERIVATIVES.  RETURNS GRADIENT IN KCAL/ANGSTROM.
C     SEE SUBROUTINE FOCK FOR MORE EXTENSIVE COMMENTS ON ENERGY TERMS.
C
C     ONLY INTERMOLECULAR TERMS ARE INCLUDED: INTRAMOLECULAR TERMS
C     ARE SET TO ZERO
C     (ASSUMED IS THAT EVERY RESIDUE CORRESPONDS TO ONE MOLECULE
C     AND VICE VERSA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      DIMENSION GXYZ(3,MAXATM)
C
C     LOCAL VARIABLES:
C
      LOGICAL ORBSI,ORBSJ,ORBSIJ,SHARE1,SHARE2
      integer resinter
      DIMENSION AGI(4),BGI(4),CGI(4),AGJ(4),BGJ(4),CGJ(4),HIJ(9,9),
     .          REPIJ(10,10),EENUC(2,10),DHII(9,9),DHJJ(9,9),DHIJ(9,9),
     .          DFII(9,9),DFJJ(9,9),DFIJ(9,9),DREPIJ(10,10),PII(9,9),
     .          PJJ(9,9),PIJ(9,9),REPIJ0(10,10),DCART(2),SCALE(2)
      DIMENSION IREPUL(4,4)
      DATA IREPUL /1,2,3,4,
     .             2,5,6,7,
     .             3,6,8,9,
     .             4,7,9,10/
      SAVE IREPUL
C
      call etimer(t1)

      DDXYZ = 1.0D-6
      DO 20 I=1,NATOMS
        DO 10 J=1,3
          GXYZ(J,I) = 0.0D0
 10     CONTINUE
 20   CONTINUE
      IF(NATOMS.EQ.1) RETURN
C
C     LOOP OVER ALL POSSIBLE PAIRS OF ATOMS.
C
      npairs = ip1(natoms+1)-1
      IJFULL = 1
      IJBOND = 1
C
C     IJFULL IS A POINTER TO POSITION IN LIST OF ALL POSSIBLE ATOM PAIRS.
C     IJBOND IS A POINTER TO POSTIION IN BONDED ATOM PAIRLIST IPAIR.
C

      do 1010 ires=2,nres
         ib = irpnt(ires)
         ie = irpnt(ires+1)-1

      DO 1000 I=ib,ie
        IAI = IATNUM(I)
        IF(IAI.EQ.0) GO TO 1000
C
C       UNPERTURBED ATOM I COORDINATES:
C
        X0 = XYZ(1,I)
        Y0 = XYZ(2,I)
        Z0 = XYZ(3,I)
C
        NORBSI = NATORB(IAI)
        ORBSI = NORBSI.GT.0
        NREPI = MAX((NORBSI*(NORBSI+1))/2,1)
        NPQI = NQUANT(IAI)
        XISI = EXPNT(0,IAI)
        XIPI = EXPNT(1,IAI)
        XIDI = EXPNT(2,IAI)
        AI0 = AL(0,IAI)
        AI1 = AL(1,IAI)
        AI2 = AL(2,IAI)
        DI1 = DL(1,IAI)
        DI2 = DL(2,IAI)
        BSI = BETA(0,IAI)
        BPI = BETA(1,IAI)
        BDI = BETA(2,IAI)
        ACI = ACORE(IAI)
        DO 40 IGAUS=1,4
          AGI(IGAUS) = AGAUS(IGAUS,IAI)
          BGI(IGAUS) = BGAUS(IGAUS,IAI)
          CGI(IGAUS) = CGAUS(IGAUS,IAI)
 40     CONTINUE
        ZNUCI = ZCHG(IAI)
C
C       TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM I.
C
        IF(ORBSI)THEN
          IJ = IIMAT(I)
          DO 70 II=1,NORBSI
            DO 60 JJ=1,II
              PII(II,JJ) = PDIAG(IJ)
              PII(JJ,II) = PDIAG(IJ)
              IJ = IJ + 1
 60         CONTINUE
 70       CONTINUE
        ENDIF
C
C       GET INTERSECTION WITH ALL OTHER ATOMS J < I.
C
        do 910 jres=1,ires-1
           jb = irpnt(jres)
           je = irpnt(jres+1)-1

        DO 900 J=jb,je
          IAJ = IATNUM(J)
          IF(IAJ.EQ.0) GO TO 900
          resinter = interres(i,j)
          NORBSJ = NATORB(IAJ)
          ORBSJ = NORBSJ.GT.0
          ORBSIJ = ORBSI.AND.ORBSJ
          call ijfind(npairs,i,j,ijaddr)
          if (ijaddr.ne.0) ijbond=ijaddr
          if(ijaddr.ne.0)then
            SHARE1 = NSHARE(1,IJBOND).NE.0
            SHARE2 = .TRUE.
          ELSE
            SHARE1 = .FALSE.
            SHARE2 = .FALSE.
          ENDIF
          NREPJ = MAX((NORBSJ*(NORBSJ+1))/2,1)
          NPQJ = NQUANT(IAJ)
          XISJ = EXPNT(0,IAJ)
          XIPJ = EXPNT(1,IAJ)
          XIDJ = EXPNT(2,IAJ)
          AJ0 = AL(0,IAJ)
          AJ1 = AL(1,IAJ)
          AJ2 = AL(2,IAJ)
          DJ1 = DL(1,IAJ)
          DJ2 = DL(2,IAJ)
          BSJ = BETA(0,IAJ)
          BPJ = BETA(1,IAJ)
          BDJ = BETA(2,IAJ)
          ACJ = ACORE(IAJ)
          DO 90 IGAUS=1,4
            AGJ(IGAUS) = AGAUS(IGAUS,IAJ)
            BGJ(IGAUS) = BGAUS(IGAUS,IAJ)
            CGJ(IGAUS) = CGAUS(IGAUS,IAJ)
 90       CONTINUE
          IF(PBC)THEN
            CALL PBCXYZ(I,J,XJ,YJ,ZJ)
          ELSE
            XJ = XYZ(1,J)
            YJ = XYZ(2,J)
            ZJ = XYZ(3,J)
          ENDIF
          ZNUCJ = ZCHG(IAJ)
C
C         TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM J.
C
          IF(ORBSJ)THEN
            IJ = IIMAT(J)
            DO 110 II=1,NORBSJ
              DO 100 JJ=1,II
                PJJ(II,JJ) = PDIAG(IJ)
                PJJ(JJ,II) = PDIAG(IJ)
                IJ = IJ + 1
 100          CONTINUE
 110        CONTINUE
          ENDIF
C
C         TEMPORARILY STORE I-J OFF-DIAGONAL BLOCK OF DENSITY MATRIX.
C
          IF(ORBSIJ.AND.SHARE1)THEN
            IJ = IJMAT(IJBOND)
            DO 130 II=1,NORBSI
              DO 120 JJ=1,NORBSJ
                PIJ(II,JJ) = PDIAT(IJ)
                IJ = IJ + 1
 120          CONTINUE
 130        CONTINUE
          ENDIF
C
C         IF THIS IS A DIRECT CALCULATION, THEN WE NEED TO COMPUTE
C         CORE-CORE AND TWO-ELECTRON INTEGRALS CORRESPONDING TO
C         THE UNPERTURBED ATOM I COORDINATES X0,Y0,Z0.  OTHERWISE
C         THEY CAN BE RETRIEVED FROM THE APPROPRIATE GLOBAL ARRAYS.
C
          IF(DIRECT.AND..NOT.CENTRL)THEN
            ILEVEL = 3
            RIJ0 = DSQRT((X0-XJ)**2 + (Y0-YJ)**2 + (Z0-ZJ)**2)
            CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,X0,Y0,Z0,ZNUCI,
     .                IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                RIJ0,ILEVEL,ENNIJ0,HIJ,REPIJ0,EENUC,resinter)
c    .                RIJ0,ILEVEL,ENNIJ0,HIJ,REPIJ0,EENUC)
          ELSEIF(.NOT.DIRECT)THEN
            IJ = IJREP(IJFULL)
            DO 220 II=1,NREPI
              DO 210 JJ=1,NREPJ
                REPIJ0(II,JJ) = EEREP(IJ)
                IJ = IJ + 1
 210          CONTINUE
 220        CONTINUE
            ENNIJ0 = ENUCLR(IJFULL)
          ENDIF
          DO 800 ICART=1,3
C
C           SHIFT EACH CARTESIAN COORDINATE OF ATOM I BY A SMALL AMOUNT
C           AND COMPUTE THE CORRESPONDING CHANGES IN DIATOMIC INTERACTIONS
C           BETWEEN ATOMS I AND J.
C
            XYZ0 = XYZ(ICART,I)
            XYZ(ICART,I) = XYZ0 + DDXYZ
            XI = XYZ(1,I)
            YI = XYZ(2,I)
            ZI = XYZ(3,I)
C
C           COMPUTE CORE-CORE, TWO-ELECTRON, AND POSSIBLY DIATOMIC
C           ONE-ELECTRON INTEGRALS FOR I-J INTERSECTION WITH PERTURBED
C           COORDINATE FOR ATOM I.  CORE-ELECTRON MAY BE SKIPPED
C           BECAUSE THEY ARE REFORMED FROM TWO-ELECTRON INTEGRALS.
C
            RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2)
            ILEVEL = 2
C
C           SKIP HIJ IF PIJ IS ZERO.
C
            IF(.NOT.SHARE1) ILEVEL = 3
C
            CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI,
     .                IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC,resinter)
c    .                RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC)
C
C           DO REVERSE STEP IF THE USER WANTS A CENTRAL DIFFERENCE.
C
            IF(CENTRL)THEN
              XYZ(ICART,I) = XYZ0 - DDXYZ
              XI = XYZ(1,I)
              YI = XYZ(2,I)
              ZI = XYZ(3,I)
              RIJOLD = RIJ
              RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2)
              CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                  DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI,
     .                  IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                  DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                  RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,EENUC,resinter)
c    .                  RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,EENUC)
C
C             COMBINE DIFFERENCE QUANTITIES FROM FORWARD AND REVERSE
C             STEPS TO GET CHANGES IN OFF-DIAGONAL BLOCK OF 1-ELECTRON
C             MATRIX AND ELECTRON REPULSIONS.
C
              DENNIJ = (ENNIJ - DENNIJ)*0.5D0
              IF(ORBSIJ.AND.SHARE1)THEN
                DO 250 II=1,NORBSI
                  DO 240 JJ=1,NORBSJ
                    DHIJ(II,JJ) = (HIJ(II,JJ) - DHIJ(II,JJ))*0.5D0
                    DFIJ(II,JJ) = DHIJ(II,JJ)
 240              CONTINUE
 250            CONTINUE
              ENDIF
C
C             NEED DREPIJ EVEN IF SPARKLES TO CALCULATE CHANGES
C             IN ELECTRON-CORE ATTRACTIONS.
C
              DO 270 II=1,NREPI
                DO 260 JJ=1,NREPJ
                  DREPIJ(II,JJ) = (REPIJ(II,JJ) - DREPIJ(II,JJ))*0.5D0
 260            CONTINUE
 270          CONTINUE
              XYZ(ICART,I) = XYZ0 + DDXYZ
              XI = XYZ(1,I)
              YI = XYZ(2,I)
              ZI = XYZ(3,I)
              RIJ = RIJOLD
            ELSE
C
C             COMPUTE CHANGES FROM FOWARD STEP ALONE.
C
              DENNIJ = ENNIJ - ENNIJ0
              IF(ORBSIJ.AND.SHARE1)THEN
                IJ = IJMAT(IJBOND)
                DO 300 II=1,NORBSI
                  DO 290 JJ=1,NORBSJ
                    DHIJ(II,JJ) = HIJ(II,JJ) - HDIAT(IJ)
                    DFIJ(II,JJ) = DHIJ(II,JJ)
                    IJ = IJ + 1
 290              CONTINUE
 300            CONTINUE
              ENDIF
              DO 320 II=1,NREPI
                DO 310 JJ=1,NREPJ
                  DREPIJ(II,JJ) = REPIJ(II,JJ) - REPIJ0(II,JJ)
 310            CONTINUE
 320          CONTINUE
            ENDIF
C
C           SINCE ELECTRON-NUCLEAR ATTRACTIONS ARE NOT STORED IN ANY
C           SORT OF GLOBAL ARRAY, THE CHANGES IN THE DIAGONAL BLOCKS
C           OF THE 1-ELECTRON MATRIX MUST BE COMPUTED DIRECTLY FROM
C           THE CHANGES IN THE ELECTRON REPULSIONS.  COPY CHANGES TO
C           FOCK MATRIX.
C
            IF(ORBSI)THEN
              DO 340 II=1,NORBSI
                DO 330 JJ=1,II
                  DHII(II,JJ) = -ZNUCJ*DREPIJ(IREPUL(II,JJ),1)
                  DFII(II,JJ) = DHII(II,JJ)
 330            CONTINUE
 340          CONTINUE
            ENDIF
            IF(ORBSJ)THEN
              DO 360 II=1,NORBSJ
                DO 350 JJ=1,II
                  DHJJ(II,JJ) = -ZNUCI*DREPIJ(1,IREPUL(II,JJ))
                  DFJJ(II,JJ) = DHJJ(II,JJ)
 350            CONTINUE
 360          CONTINUE
            ENDIF
C
C           ADD ON 2-CENTER 2-ELECTRON TERMS TO THE CHANGES IN THE
C           OFF-DIAGONAL BLOCK OF THE FOCK MATRIX
C
            IF(ORBSIJ.AND.SHARE1)THEN
              DO 400 IG=1,NORBSI
                DO 390 JG=1,NORBSJ
                  DGIJ = 0.0D0
                  DO 380 IP=1,NORBSI
                    IR = IREPUL(IG,IP)
                    DO 370 JP=1,NORBSJ
                      JR = IREPUL(JG,JP)
                      DGIJ = DGIJ + PIJ(IP,JP)*DREPIJ(IR,JR)
 370                CONTINUE
 380              CONTINUE
                  DFIJ(IG,JG) = DFIJ(IG,JG) - 0.5D0*DGIJ
 390            CONTINUE
 400          CONTINUE
            ENDIF
C
C           ADD ON 2-CENTER 2-ELECTRON TERMS TO CHANGES IN DIAGONAL
C           BLOCKS OF FOCK MATRIX.
C
            IF(ORBSIJ)THEN
              DO 440 IG=1,NORBSI
                DO 430 JG=1,IG
                  JJG = IREPUL(JG,JG)
                  IJG = IREPUL(IG,JG)
                  DGIJ = 0.0D0
                  DO 420 IP=1,NORBSJ
                    DO 410 JP=1,NORBSJ
                      IJP = IREPUL(IP,JP)
                      DGIJ = DGIJ + DREPIJ(IJG,IJP)*PJJ(IP,JP)
 410                CONTINUE
 420              CONTINUE
                  DFII(IG,JG) = DFII(IG,JG) + DGIJ
 430            CONTINUE
 440          CONTINUE
              DO 480 IG=1,NORBSJ
                DO 470 JG=1,IG
                  JJG = IREPUL(JG,JG)
                  IJG = IREPUL(IG,JG)
                  DGIJ = 0.0D0
                  DO 460 IP=1,NORBSI
                    DO 450 JP=1,NORBSI
                      IJP = IREPUL(IP,JP)
                      DGIJ = DGIJ + DREPIJ(IJP,IJG)*PII(IP,JP)
 450                CONTINUE
 460              CONTINUE
                  DFJJ(IG,JG) = DFJJ(IG,JG) + DGIJ
 470            CONTINUE
 480          CONTINUE
            ENDIF
C
C           WE NOW HAVE ALL CHANGES IN THE 1-ELECTRON AND FOCK MATRICES
C           DUE TO THE I-J INTERACTION.  NOW COMPUTE THE CORRESPONDING
C           CHANGE IN THE TOTAL ENERGY.
C
            DIJ = 0.0D0
            DII = 0.0D0
            DJJ = 0.0D0
            IF(ORBSI)THEN
              DO 520 II=1,NORBSI
                DO 510 JJ=1,II
                  DDIJ = (DHII(II,JJ) + DFII(II,JJ))*PII(II,JJ)
                  DIJ = DIJ + DDIJ
 510            CONTINUE
                DII = DII + DDIJ
 520          CONTINUE
            ENDIF
            IF(ORBSJ)THEN
              DO 540 II=1,NORBSJ
                DO 530 JJ=1,II
                  DDIJ = (DHJJ(II,JJ) + DFJJ(II,JJ))*PJJ(II,JJ)
                  DIJ = DIJ + DDIJ
 530            CONTINUE
                DJJ = DJJ + DDIJ
 540          CONTINUE
            ENDIF
            DELECT = DIJ - 0.5D0*(DII + DJJ)
            IF(ORBSIJ.AND.SHARE1)THEN
              DO 560 II=1,NORBSI
                DO 550 JJ=1,NORBSJ
                  DELECT = DELECT + (DHIJ(II,JJ)+DFIJ(II,JJ))*PIJ(II,JJ)
 550            CONTINUE
 560          CONTINUE
            ENDIF
C
C           UPDATE GRADIENT ENTRIES FOR ATOMS I AND J:
C
            DGRAD1 = DELECT + DENNIJ
            GXYZ(ICART,I) = GXYZ(ICART,I) + DGRAD1
            GXYZ(ICART,J) = GXYZ(ICART,J) - DGRAD1
C
C           RESTORE CARTESIAN COORDINATE OF ATOM I.
C
            XYZ(ICART,I) = XYZ0
 800      CONTINUE
C
C         UPDATE POINTERS FOR FULL AND BOND-BASED PAIRLIST.
C
          IJFULL = IJFULL + 1
 900    CONTINUE
 910  continue
 1000 CONTINUE
 1010 continue

#if 0
C     the following code (between if 0, endif directives) will turn
C     the code into one that will include both intramolecular terms
C     and intermolecular terms

      IJFULL = 1
      IJBOND = 1
C
C     IJFULL IS A POINTER TO POSITION IN LIST OF ALL POSSIBLE ATOM PAIRS.
C     IJBOND IS A POINTER TO POSTIION IN BONDED ATOM PAIRLIST IPAIR.
C
      do 2010 ires=1,nres
         ib = irpnt(ires)
         ie = irpnt(ires+1)-1

      DO 2000 I=ib+1,ie
        IAI = IATNUM(I)
        IF(IAI.EQ.0) GO TO 2000
C
C       UNPERTURBED ATOM I COORDINATES:
C
        X0 = XYZ(1,I)
        Y0 = XYZ(2,I)
        Z0 = XYZ(3,I)
C
        NORBSI = NATORB(IAI)
        ORBSI = NORBSI.GT.0
        NREPI = MAX((NORBSI*(NORBSI+1))/2,1)
        NPQI = NQUANT(IAI)
        XISI = EXPNT(0,IAI)
        XIPI = EXPNT(1,IAI)
        XIDI = EXPNT(2,IAI)
        AI0 = AL(0,IAI)
        AI1 = AL(1,IAI)
        AI2 = AL(2,IAI)
        DI1 = DL(1,IAI)
        DI2 = DL(2,IAI)
        BSI = BETA(0,IAI)
        BPI = BETA(1,IAI)
        BDI = BETA(2,IAI)
        ACI = ACORE(IAI)
        DO 1040 IGAUS=1,4
          AGI(IGAUS) = AGAUS(IGAUS,IAI)
          BGI(IGAUS) = BGAUS(IGAUS,IAI)
          CGI(IGAUS) = CGAUS(IGAUS,IAI)
 1040  CONTINUE
        ZNUCI = ZCHG(IAI)
C
C       TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM I.
C
        IF(ORBSI)THEN
          IJ = IIMAT(I)
          DO 1070 II=1,NORBSI
            DO 1060 JJ=1,II
              PII(II,JJ) = PDIAG(IJ)
              PII(JJ,II) = PDIAG(IJ)
              IJ = IJ + 1
 1060         CONTINUE
 1070       CONTINUE
        ENDIF
C
C       GET INTERSECTION WITH ALL OTHER ATOMS J < I.
C

        DO 1900 J=ib,i-1
          IAJ = IATNUM(J)
          IF(IAJ.EQ.0) GO TO 1900
          resinter = interres(i,j)
          NORBSJ = NATORB(IAJ)
          ORBSJ = NORBSJ.GT.0
          ORBSIJ = ORBSI.AND.ORBSJ
          call ijfind(npairs,i,j,ijaddr)
          if (ijaddr.ne.0) ijbond=ijaddr
          if(ijaddr.ne.0)then
            SHARE1 = NSHARE(1,IJBOND).NE.0
            SHARE2 = .TRUE.
          ELSE
            SHARE1 = .FALSE.
            SHARE2 = .FALSE.
          ENDIF
          NREPJ = MAX((NORBSJ*(NORBSJ+1))/2,1)
          NPQJ = NQUANT(IAJ)
          XISJ = EXPNT(0,IAJ)
          XIPJ = EXPNT(1,IAJ)
          XIDJ = EXPNT(2,IAJ)
          AJ0 = AL(0,IAJ)
          AJ1 = AL(1,IAJ)
          AJ2 = AL(2,IAJ)
          DJ1 = DL(1,IAJ)
          DJ2 = DL(2,IAJ)
          BSJ = BETA(0,IAJ)
          BPJ = BETA(1,IAJ)
          BDJ = BETA(2,IAJ)
          ACJ = ACORE(IAJ)
          DO 1090 IGAUS=1,4
            AGJ(IGAUS) = AGAUS(IGAUS,IAJ)
            BGJ(IGAUS) = BGAUS(IGAUS,IAJ)
            CGJ(IGAUS) = CGAUS(IGAUS,IAJ)
 1090       CONTINUE
          IF(PBC)THEN
            CALL PBCXYZ(I,J,XJ,YJ,ZJ)
          ELSE
            XJ = XYZ(1,J)
            YJ = XYZ(2,J)
            ZJ = XYZ(3,J)
          ENDIF
          ZNUCJ = ZCHG(IAJ)
C
C         TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM J.
C
          IF(ORBSJ)THEN
            IJ = IIMAT(J)
            DO 1110 II=1,NORBSJ
              DO 1100 JJ=1,II
                PJJ(II,JJ) = PDIAG(IJ)
                PJJ(JJ,II) = PDIAG(IJ)
                IJ = IJ + 1
 1100          CONTINUE
 1110        CONTINUE
          ENDIF
C
C         TEMPORARILY STORE I-J OFF-DIAGONAL BLOCK OF DENSITY MATRIX.
C
          IF(ORBSIJ.AND.SHARE1)THEN
            IJ = IJMAT(IJBOND)
            DO 1130 II=1,NORBSI
              DO 1120 JJ=1,NORBSJ
                PIJ(II,JJ) = PDIAT(IJ)
                IJ = IJ + 1
 1120          CONTINUE
 1130        CONTINUE
          ENDIF
C
C         IF THIS IS A DIRECT CALCULATION, THEN WE NEED TO COMPUTE
C         CORE-CORE AND TWO-ELECTRON INTEGRALS CORRESPONDING TO
C         THE UNPERTURBED ATOM I COORDINATES X0,Y0,Z0.  OTHERWISE
C         THEY CAN BE RETRIEVED FROM THE APPROPRIATE GLOBAL ARRAYS.
C
          IF(DIRECT.AND..NOT.CENTRL)THEN
            ILEVEL = 3
            RIJ0 = DSQRT((X0-XJ)**2 + (Y0-YJ)**2 + (Z0-ZJ)**2)
            CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,X0,Y0,Z0,ZNUCI,
     .                IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                RIJ0,ILEVEL,ENNIJ0,HIJ,REPIJ0,EENUC,resinter)
c    .                RIJ0,ILEVEL,ENNIJ0,HIJ,REPIJ0,EENUC)
          ELSEIF(.NOT.DIRECT)THEN
            IJ = IJREP(IJFULL)
            DO 1220 II=1,NREPI
              DO 1210 JJ=1,NREPJ
                REPIJ0(II,JJ) = EEREP(IJ)
                IJ = IJ + 1
 1210          CONTINUE
 1220        CONTINUE
            ENNIJ0 = ENUCLR(IJFULL)
          ENDIF
          DO 1800 ICART=1,3
C
C           SHIFT EACH CARTESIAN COORDINATE OF ATOM I BY A SMALL AMOUNT
C           AND COMPUTE THE CORRESPONDING CHANGES IN DIATOMIC INTERACTIONS
C           BETWEEN ATOMS I AND J.
C
            XYZ0 = XYZ(ICART,I)
            XYZ(ICART,I) = XYZ0 + DDXYZ
            XI = XYZ(1,I)
            YI = XYZ(2,I)
            ZI = XYZ(3,I)
C
C           COMPUTE CORE-CORE, TWO-ELECTRON, AND POSSIBLY DIATOMIC
C           ONE-ELECTRON INTEGRALS FOR I-J INTERSECTION WITH PERTURBED
C           COORDINATE FOR ATOM I.  CORE-ELECTRON MAY BE SKIPPED
C           BECAUSE THEY ARE REFORMED FROM TWO-ELECTRON INTEGRALS.
C
            RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2)
            ILEVEL = 2
C
C           SKIP HIJ IF PIJ IS ZERO.
C
            IF(.NOT.SHARE1) ILEVEL = 3
C
            CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI,
     .                IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC,resinter)
c    .                RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC)
C
C           DO REVERSE STEP IF THE USER WANTS A CENTRAL DIFFERENCE.
C
            IF(CENTRL)THEN
              XYZ(ICART,I) = XYZ0 - DDXYZ
              XI = XYZ(1,I)
              YI = XYZ(2,I)
              ZI = XYZ(3,I)
              RIJOLD = RIJ
              RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2)
              CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                  DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI,
     .                  IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                  DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .                  RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,EENUC,resinter)
c    .                  RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,EENUC)
C
C             COMBINE DIFFERENCE QUANTITIES FROM FORWARD AND REVERSE
C             STEPS TO GET CHANGES IN OFF-DIAGONAL BLOCK OF 1-ELECTRON
C             MATRIX AND ELECTRON REPULSIONS.
C
              DENNIJ = (ENNIJ - DENNIJ)*0.5D0
              IF(ORBSIJ.AND.SHARE1)THEN
                DO 1250 II=1,NORBSI
                  DO 1240 JJ=1,NORBSJ
                    DHIJ(II,JJ) = (HIJ(II,JJ) - DHIJ(II,JJ))*0.5D0
                    DFIJ(II,JJ) = DHIJ(II,JJ)
 1240              CONTINUE
 1250            CONTINUE
              ENDIF
C
C             NEED DREPIJ EVEN IF SPARKLES TO CALCULATE CHANGES
C             IN ELECTRON-CORE ATTRACTIONS.
C
              DO 1270 II=1,NREPI
                DO 1260 JJ=1,NREPJ
                  DREPIJ(II,JJ) = (REPIJ(II,JJ) - DREPIJ(II,JJ))*0.5D0
 1260            CONTINUE
 1270          CONTINUE
              XYZ(ICART,I) = XYZ0 + DDXYZ
              XI = XYZ(1,I)
              YI = XYZ(2,I)
              ZI = XYZ(3,I)
              RIJ = RIJOLD
            ELSE
C
C             COMPUTE CHANGES FROM FOWARD STEP ALONE.
C
              DENNIJ = ENNIJ - ENNIJ0
              IF(ORBSIJ.AND.SHARE1)THEN
                IJ = IJMAT(IJBOND)
                DO 1300 II=1,NORBSI
                  DO 1290 JJ=1,NORBSJ
                    DHIJ(II,JJ) = HIJ(II,JJ) - HDIAT(IJ)
                    DFIJ(II,JJ) = DHIJ(II,JJ)
                    IJ = IJ + 1
 1290              CONTINUE
 1300            CONTINUE
              ENDIF
              DO 1320 II=1,NREPI
                DO 1310 JJ=1,NREPJ
                  DREPIJ(II,JJ) = REPIJ(II,JJ) - REPIJ0(II,JJ)
 1310            CONTINUE
 1320          CONTINUE
            ENDIF
C
C           SINCE ELECTRON-NUCLEAR ATTRACTIONS ARE NOT STORED IN ANY
C           SORT OF GLOBAL ARRAY, THE CHANGES IN THE DIAGONAL BLOCKS
C           OF THE 1-ELECTRON MATRIX MUST BE COMPUTED DIRECTLY FROM
C           THE CHANGES IN THE ELECTRON REPULSIONS.  COPY CHANGES TO
C           FOCK MATRIX.
C
            IF(ORBSI)THEN
              DO 1340 II=1,NORBSI
                DO 1330 JJ=1,II
                  DHII(II,JJ) = -ZNUCJ*DREPIJ(IREPUL(II,JJ),1)
                  DFII(II,JJ) = DHII(II,JJ)
 1330            CONTINUE
 1340          CONTINUE
            ENDIF
            IF(ORBSJ)THEN
              DO 1360 II=1,NORBSJ
                DO 1350 JJ=1,II
                  DHJJ(II,JJ) = -ZNUCI*DREPIJ(1,IREPUL(II,JJ))
                  DFJJ(II,JJ) = DHJJ(II,JJ)
 1350            CONTINUE
 1360          CONTINUE
            ENDIF
C
C           ADD ON 2-CENTER 2-ELECTRON TERMS TO THE CHANGES IN THE
C           OFF-DIAGONAL BLOCK OF THE FOCK MATRIX
C
            IF(ORBSIJ.AND.SHARE1)THEN
              DO 1400 IG=1,NORBSI
                DO 1390 JG=1,NORBSJ
                  DGIJ = 0.0D0
                  DO 1380 IP=1,NORBSI
                    IR = IREPUL(IG,IP)
                    DO 1370 JP=1,NORBSJ
                      JR = IREPUL(JG,JP)
                      DGIJ = DGIJ + PIJ(IP,JP)*DREPIJ(IR,JR)
 1370                CONTINUE
 1380              CONTINUE
                  DFIJ(IG,JG) = DFIJ(IG,JG) - 0.5D0*DGIJ
 1390            CONTINUE
 1400          CONTINUE
            ENDIF
C
C           ADD ON 2-CENTER 2-ELECTRON TERMS TO CHANGES IN DIAGONAL
C           BLOCKS OF FOCK MATRIX.
C
            IF(ORBSIJ)THEN
              DO 1440 IG=1,NORBSI
                DO 1430 JG=1,IG
                  JJG = IREPUL(JG,JG)
                  IJG = IREPUL(IG,JG)
                  DGIJ = 0.0D0
                  DO 1420 IP=1,NORBSJ
                    DO 1410 JP=1,NORBSJ
                      IJP = IREPUL(IP,JP)
                      DGIJ = DGIJ + DREPIJ(IJG,IJP)*PJJ(IP,JP)
 1410                CONTINUE
 1420              CONTINUE
                  DFII(IG,JG) = DFII(IG,JG) + DGIJ
 1430            CONTINUE
 1440          CONTINUE
              DO 1480 IG=1,NORBSJ
                DO 1470 JG=1,IG
                  JJG = IREPUL(JG,JG)
                  IJG = IREPUL(IG,JG)
                  DGIJ = 0.0D0
                  DO 1460 IP=1,NORBSI
                    DO 1450 JP=1,NORBSI
                      IJP = IREPUL(IP,JP)
                      DGIJ = DGIJ + DREPIJ(IJP,IJG)*PII(IP,JP)
 1450                CONTINUE
 1460              CONTINUE
                  DFJJ(IG,JG) = DFJJ(IG,JG) + DGIJ
 1470            CONTINUE
 1480          CONTINUE
            ENDIF
C
C           WE NOW HAVE ALL CHANGES IN THE 1-ELECTRON AND FOCK MATRICES
C           DUE TO THE I-J INTERACTION.  NOW COMPUTE THE CORRESPONDING
C           CHANGE IN THE TOTAL ENERGY.
C
            DIJ = 0.0D0
            DII = 0.0D0
            DJJ = 0.0D0
            IF(ORBSI)THEN
              DO 1520 II=1,NORBSI
                DO 1510 JJ=1,II
                  DDIJ = (DHII(II,JJ) + DFII(II,JJ))*PII(II,JJ)
                  DIJ = DIJ + DDIJ
 1510            CONTINUE
                DII = DII + DDIJ
 1520          CONTINUE
            ENDIF
            IF(ORBSJ)THEN
             DO 1540 II=1,NORBSJ
                DO 1530 JJ=1,II
                  DDIJ = (DHJJ(II,JJ) + DFJJ(II,JJ))*PJJ(II,JJ)
                  DIJ = DIJ + DDIJ
 1530            CONTINUE
                DJJ = DJJ + DDIJ
 1540          CONTINUE
            ENDIF
            DELECT = DIJ - 0.5D0*(DII + DJJ)
            IF(ORBSIJ.AND.SHARE1)THEN
              DO 1560 II=1,NORBSI
                DO 1550 JJ=1,NORBSJ
                  DELECT = DELECT + (DHIJ(II,JJ)+DFIJ(II,JJ))*PIJ(II,JJ)
 1550            CONTINUE
 1560          CONTINUE
            ENDIF
C
C           UPDATE GRADIENT ENTRIES FOR ATOMS I AND J:
C
            DGRAD1 = DELECT + DENNIJ
            GXYZ(ICART,I) = GXYZ(ICART,I) + DGRAD1
            GXYZ(ICART,J) = GXYZ(ICART,J) - DGRAD1
C
C           RESTORE CARTESIAN COORDINATE OF ATOM I.
C
            XYZ(ICART,I) = XYZ0
 1800      CONTINUE
C
C         UPDATE POINTERS FOR FULL AND BOND-BASED PAIRLIST.
C
          IJFULL = IJFULL + 1
 1900    CONTINUE
 2000 CONTINUE
 2010 continue
#endif

C
C     CONVERT ENERGY DIFFERENCES TO DERIVATIVES IN KCAL/ANGSTROM.
C
      FACTR = 23.061D0/DDXYZ
      DO 2400 I=1,NATOMS
        DO 2300 J=1,3
          GXYZ(J,I) = GXYZ(J,I)*FACTR
 2300   CONTINUE
 2400 CONTINUE

      call etimer(t2)
      tgrad = t2-t1
      RETURN
      END