C  4 MAR 08 - MWS - ADD DOWNSHIFT TO C1 SYMMETRY
C 20 AUG 07 - KI  - NEW PARALLEL PROGRAM FROM IMS
C
C*MODULE MP2GR2  *DECK MP2GRDIMS
      SUBROUTINE MP2GRDIMS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /FMCOM/ X(1)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      NPROCS = NPROC
      NSHLLS = NSHELL
C
      CALL VALFM(LOADFM)
      LNUMP  = LOADFM + 1
      LNLOOP = LNUMP  + NPROCS + 1
      LILOOP = LNLOOP + NPROCS*4
      LLENDD = LILOOP + NSHLLS*NPROCS*2
      LLRECV = LLENDD + NPROCS*3
      LAST   = LLRECV + NPROCS*2
      NEED = LAST - LOADFM - 1
C
      CALL GETFM(NEED)
      CALL MP2GRDIMSX(X(LNUMP),X(LNLOOP),X(LILOOP),X(LLENDD),
     *                X(LLRECV),NPROCS,NSHLLS)
      CALL RETFM(NEED)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MP2GRDIMSX
      SUBROUTINE MP2GRDIMSX(NUMP,NLOOP,ILOOP,LENDDI,LRECV,NPROCS,NSHLLS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DSKSAV,NXT,DIVANG
C
      PARAMETER (MXSH=5000, MXGTOT=20000,
     *           MXRT=100, MXATM=2000)
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
C
      COMMON /ENRGMP/ EMP2,EMP3,EMP4
      COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /INFOA / NAT,ICH,MUL,MNUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
C
      DIMENSION NUMP(0:NPROCS),NLOOP(NPROCS,4),ILOOP(NSHLLS,NPROCS,2),
     *          LENDDI(0:NPROCS-1,3),LRECV(0:NPROCS-1,2)
C
      DATA CHECK/8HCHECK   /
C     DATA NFT9,NFT15,NFT16,NFT19,NFT20/9,15,16,19,20/
      DATA NFT9,NFT15,NFT16,NFT19,NFT20/50,51,52,53,54/
C
C     ----- IMS CLOSED-SHELL MP2 ENERGY GRADIENT DRIVER -----
C
      QQ4=ONE
      NXT = IBTYP.EQ.1
      E2= ZER
      ITERMX= MAXITC
      CUTOFF=TOL
C
      DSKSAV=DSKWRK
      DSKWRK=.TRUE.
C
C     SUBTRACT FMO PROJECTED OUT ORBITALS.
C     THIS DOES NOT WORK. IT IS PROBABLY RELATED TO THE CPHF EQUATIONS
C     HAVING TO BE SOLVED IN THE FULL SPACE (OR IS IT A PROGRAM
C     LIMITATION THAT NVIR MUST BE EQUAL TO NVAC?)
C     NOTE THAT NOT SUBTRACTING NORBPROJ LEADS TO A TRIFLE DIFFERENCE IN
C     THE RESULTS (1E-6 A.U. OR SMALLER) AND SOME TIME INCREASE.
C     IF(NFG.NE.0) NO = NO - NORBPROJ
      CALL SETFMPI(MEMPRI)
C
C    ----- SIZE OF DISK CACHE -----
C
      NCACHE= 1000000
C
C    ----- VARIABLES FOR MP2 ENERGY GRADIENT -----
C
C NO(=NQMT): NUMBER OF CANONICAL ORBITALS
C NBF      : NUMBER OF BASIS FUNCTIONS
C NSHELL   : NUMBER OF SHELLS
C NOCC     : NUMBER OF OCCUPIED ORBITALS
C NACORE   : NUMBER OF CORE ORBITALS
C NOAC     : NUMBER OF ACTIVE (CORRELATED) ORBITALS
C NVIR     : NUMBER OF VIRTUAL ORBITALS
C NVAC     : NUMBER OF ACTIVE VIRTUAL ORBITALS
C NOC2     : NUMBER OF IJ PAIR
C
      NOCC = NOA
      NVIR = NQMT  - NOA
      NOAC = NOA - NACORE
      NVAC = NQMT  - NOA
C     NVAC = NO  - NOA
C
      NBF2 = NBF*NBF
      NBF3 =(NBF*NBF+NBF)/2
      NSH2 =(NSHELL*NSHELL+NSHELL)/2
      NOCC2= NOCC*NOCC
      NOAC2= NOAC*NOAC
      NVIR2= NVIR*NVIR
      NOCC3=(NOCC*NOCC+NOCC)/2
C
      IF(MASWRK) WRITE(IW,9000) NOCC,NVIR,NOAC,NVAC,NBF,NSHELL,CUTOFF
C     IF(MASWRK.AND.NFG.NE.0.AND.NORBPROJ.NE.0) WRITE(IW,9005) NORBPROJ
C
C          THIS CODE WORKS ONLY IN THE C1 POINT GROUP AT PRESENT.
      IF(NT.GT.1  .AND.  MASWRK) WRITE(IW,9004)
 9004 FORMAT(/1X,'THE USE OF POINT GROUP SYMMETRY HAS BEEN TURNED OFF.')
      CALL SYMOFF
C
      CALL BASCHK(LMAX)
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      IF(LMAX.EQ.5) NANGM = 21
      IF(LMAX.EQ.6) NANGM = 28
      MAXG = NANGM**4
C
C        RUN THROUGH THE PROGRAM AS A SIMULATED 'CHECK' RUN,
C        THEN COME BACK TO THIS POINT IF IT IS A REAL RUN.
C
      EXESAV = EXETYP
      EXETYP = CHECK
      ICHKPASS=0
   50 CONTINUE
      ICHKPASS=ICHKPASS+1
C
C     ----- PRINT OUT NEEDED MEMORY AND DISK SPACE -----
C     THE MEMORY REQUIREMENTS JUST BELOW DON'T MATCH THE ACTUAL CODE.
C
C---      NEEDA  = NBF2+3*NOCC*NBF+3*NVIR*NBF+NBF+MAXG+NSH2
C---      NEEDBMX= (NANGM**3)*(NBF+NOCC)+NCACHE+NANGM*NOCC2*NBF
C---      NEEDBMN= (NANGM**3)*(NBF+NOCC)+NCACHE+      NOCC2*NBF
C---      NEEDC  = 2*(NOCC2+NVIR2+NBF2)+3*NOCC*NVIR
C---      NEEDD  = NOCC2*NBF+NOCC*NOAC+NOAC2+NBF*(MAX(NOCC,NVIR))
C---      NEEDA2 = 2*NOCC2*NBF
C---      NEEDB2 = 2*NOAC2*(NVAC+NBF)
C---      NEEDEMX= (NANGM**3)*NBF+(NANGM**2)*NOCC*NBF
C---     *        +NOAC2*NBF*(1+NANGM)
C---     *        +2*NOCC*NBF+NCACHE+NANGM*NOCC3*NBF
C---      NEEDEMN= (NANGM**3)*NBF+(NANGM**2)*NOCC*NBF
C---     *        +NOAC2*NBF*(1+NANGM)
C---     *        +2*NOCC*NBF+NCACHE+      NOCC3*NBF
C---      NEEDF  = 2*NOCC3*(NBF+NVIR)+NOCC3*NOCC+NOCC3+NOCC2+NOCC*NVIR
C---     *        +NCACHE
C---      NEED1 = NEEDA+NEEDBMX
C---      NEED2 = NEEDA+NEEDBMN
C---      NEED3 = NEEDA+NEEDC+NEEDD+NEEDA2
C---      NEED4 = NEEDA+NEEDC+NEEDD+NEEDB2
C---      NEED5 = NEEDA+NEEDC+NEEDEMX
C---      NEED6 = NEEDA+NEEDC+NEEDEMN
C---      NEED7 = NEEDA+NEEDC+NEEDF
C---      MXNEED= MAX(NEED1,NEED3)
C---      MXNEED= MAX(MXNEED,NEED4)
C---      MXNEED= MAX(MXNEED,NEED5)
C---      MXNEED= MAX(MXNEED,NEED7)
C---      MINEED= MAX(NEED2,NEED3)
C---      MINEED= MAX(MINEED,NEED4)
C---      MINEED= MAX(MINEED,NEED6)
C---      MINEED= MAX(MINEED,NEED7)
C---C
C---      IF(MASWRK) WRITE(IW,9010) MXNEED,(MXNEED-1)/1000000+1,
C---     *                          MINEED,(MINEED-1)/1000000+1
C
      RNOCC = NOCC
      RNVIR = NVIR
      RNBF  = NBF
      RNOAC = NOAC
      RNVAC = NVAC
      DNEED1= RNOCC*RNOCC*RNVIR*(RNBF+RNVIR+RNOCC)+
     *        RNOAC*RNOAC*RNVAC*RNBF
      DNEED2= RNOAC*RNOAC*RNBF*RNBF+RNOCC*RNOCC*RNBF*(RNVIR+RNBF/2.0D+0)
      MXDISK= INT(MAX(DNEED1,DNEED2)/125000.0D+0)
      IF(GOPARR) MXDISK= MXDISK/NPROC
      IF(MASWRK  .AND.  ICHKPASS.EQ.1) WRITE(IW,9015) MXDISK
C
C---      CALL GOTFMP(NREQ)
C---      IF(NREQ.LT.MINEED) THEN
C---         WRITE(IW,9020) NREQ/1000000
C---         CALL ABRT
C---      ENDIF
C
C     ----- DISTRIBUTE P OF (PQ|RS) -----
C     ----- FOR STATIC LOAD BALANCING -----
C
      IF((.NOT.NXT).AND.GOPARR) THEN
         IPROC= 1
         IFUNC= NBF/NPROC
         NUMP(0)= 0
         DO I=1,NPROC
            NUMP(I)= NSHELL
            NLOOP(I,2)= 0
         ENDDO
         DO I=1,NSHELL
            IF(KLOC(I)-1.GE.IFUNC) THEN
               NUMP(IPROC)= I-1
               IPROC= IPROC+1
               IFUNC= (NBF*IPROC)/NPROC
            ENDIF
         ENDDO
         IF(NUMP(IPROC-1).NE.NSHELL) THEN
            MXP=NBF-KLOC(NUMP(IPROC-1)+1)+1
            NLOOP(IPROC,2)=MXP
         ENDIF
         DO I=1,IPROC-1
            IF(NUMP(I).NE.NSHELL) THEN
               IMXP= KLOC(NUMP(I)+1)-KLOC(NUMP(I-1)+1)
               NLOOP(I,2)=IMXP
               IF(IMXP.GT.MXP) MXP= IMXP
            ENDIF
         ENDDO
C
C            DYNAMIC LOAD BALANCING SETS REAL NLOOP VALUES IN -MPGQTA-
C            SERIAL DOESN'T SEEM TO NEED TO HAVE NLOOP SET TO ANYTHING
C
      ELSE
         DO IPROC=1,NPROCS
            DO I=1,4
               NLOOP(IPROC,I) = 0
            ENDDO
         ENDDO
      ENDIF
C
C     ----- DISTRIBUTE P OF (PQ|RS) -----
C     ----- FOR BACK-TRANSFORMATION -----
C
      IF(GOPARR) THEN
         IPROC=0
         NEXT=1
         ICOUNT=1
         DO II= 1,NPROC
            NLOOP(II,3)= 0
            NLOOP(II,4)= 0
         ENDDO
         DO ISH= NSHELL,1,-1
            IPROC= IPROC+NEXT
            IF(IPROC.GT.NPROC) THEN
               IPROC= NPROC
               NEXT =-1
               ICOUNT=ICOUNT+1
            ENDIF
            IF(IPROC.LT.1) THEN
               IPROC= 1
               NEXT = 1
               ICOUNT= ICOUNT+1
            ENDIF
            ILOOP(ICOUNT,IPROC,2)= ISH
            NLOOP(IPROC,3)=NLOOP(IPROC,3)+1
            NLOOP(IPROC,4)=NLOOP(IPROC,4)+KMAX(ISH)-KMIN(ISH)+1
         ENDDO
      ENDIF
C
C    ----- ALLOCATE MEMORY FOR MP2 ENERGY GRADIENT -----
C
C IVEC    : MATRIX OF MO COEFFICIENTS        NBF*NBF
C IVECO   : MATRIX OF OCCUPIED               NOCC*NBF
C           MO COEEFICIENTS
C IVECV   : MATRIX OF VIRTUAL                NVIR*NBF
C           MO COEFFICIENTS
C IENG    : ORBITAL ENERGIES                 NBF
C IGHND   : AO INTEGRALS FOR D,F,G           MAXG --> NANGM**4
C IXINT   : EXCHANGE INTEGRALS               (NSHELL**2 + NSHELL)/2
C
      CALL VALFMP(LOADFM)
      IVEC    = 1       + LOADFM
      IVECO   = IVEC    + NBF2
      IVECO2  = IVECO   + NOCC*NBF
      IVECO3  = IVECO2  + NOCC*NBF
      IVECV   = IVECO3  + NOCC*NBF
      IVECV2  = IVECV   + NVIR*NBF
      IVECV3  = IVECV2  + NVIR*NBF
      IENG    = IVECV3  + NVIR*NBF
      IGHND   = IENG    + NBF
      IXINT   = IGHND   + MAXG
      ILBLINT = IXINT   + NSH2
      INUMINT = ILBLINT + NBF*NANGM**3
      LAST    = INUMINT + NANGM**3
      NEEDA = LAST  - LOADFM -1
      CALL GETFM(NEEDA)
C
C     ----- READ MO COEFFICIENTS AND MO ENERGIES ------
C
      CALL DAREAD(IDAF,IODA,X(IVEC),NBF2,15,0)
      CALL DAREAD(IDAF,IODA,X(IENG),NBF ,17,0)
C
      MAXC = IDAMAX(NBF*NQMT,X(IVEC),1)
      IF(ABS(X(IVEC+MAXC-1)).GT.1.0D+02) THEN
         MAXM = 1 + MAXC/NBF
         MAXA = MAXC - NBF*(MAXM-1)
         IF(MASWRK) WRITE(IW,9030) ABS(X(IVEC+MAXC-1)),MAXA,MAXM
      END IF
C
C     ----- READ SCHWARTZ INEQUALITY INFORMATION -----
C
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(X(IXINT),1,NSH2)
      ELSE
         IF(ISCHWZ.EQ.1) THEN
            CALL DAREAD(IDAF,IODA,X(IXINT),NSH2,54,0)
         ELSE
            CALL DCOPY(NSH2,1.0D+01,0,X(IXINT),1)
         END IF
      END IF
C
C     ----- ALLOCATE MEMORY FOR 1ST, 2ND, 3RD TRANSFORMATION -----
C
C IINT    : AO INTEGRALS                     (NANGM**3)*NBF
C IPTR1   : PARTIALLY TRANSFORMED INTEGRAL   (NANGM**3)*NOC
C IPTR2   : PARTIALLY TRANSFORMED INTEGRAL   NANGM*NOCC*NOCC*NBF
C IPTR3   : PARTIALLY TRANSFORMED INTEGRAL   NOCC*NOCC*NB3
C NB3     : NUMBER OF VIRTUAL ORBITAL B IN A BATCH AT 3RD TRANS.
C
      NB3= NCACHE/NOCC2 + 1
      IF(NB3.GT.NVIR) NB3= NVIR
C
      CALL GOTFMP(NPTR2)
      NPTR2= NPTR2-(NANGM**3)*NBF-(NANGM**3)*NOCC-NOCC2*NB3
      JPTR2 = NANGM*NOCC2*NBF
      DIVANG=.FALSE.
      IF(NPTR2.LT.JPTR2) THEN
         DO 110 I=1, NANGM-1
            MXANGM= NANGM-I
            JPTR2= MXANGM*NOCC2*NBF
            IF(NPTR2.GE.JPTR2) GO TO 120
  110    CONTINUE
         IF(MASWRK) WRITE(IW,9040) (NEEDA+JPTR2)/1000000
         CALL ABRT
  120    DIVANG=.TRUE.
         IF(MASWRK) WRITE(IW,9050) MXANGM
      ENDIF
C
      CALL VALFMP(LOADFM)
      IINT  = 1     + LOADFM
      IPTR1 = IINT  +(NANGM**3)*NBF
      IPTR2 = IPTR1 +(NANGM**3)*NOCC
      IPTR3 = IPTR2 + JPTR2
      LAST  = IPTR3 + NOCC2*NB3
      NEEDB = LAST  - LOADFM -1
      CALL GETFMP(NEEDB)
      NEED2 = NEEDA + NEEDB
C
      IF(EXETYP.EQ.CHECK) GO TO 130
C
C     ----- OPEN DIRECT ACCESS FILE (WORK15) -----
C
      CALL MPOPEN(NFT15,NOCC2,'MCQD51')
C
C     ----- 1ST, 2ND, 3RD QUARTER TRANSFORMATIONS  -----
C
      CALL MOCOPY(X(IVEC),X(IVECO),NOCC,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,3)
      CALL MOCOPY(X(IVEC),X(IVECV),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,2)
      CALL MPGQTA(NFT15,X(IVECO),X(IVECV),X(IGHND),X(IXINT),
     *            X(IINT),X(IPTR1),X(IPTR2),X(IPTR3),
     *            NOCC,NVIR,NBF,NOAC,NANGM,CUTOFF,NUMP,
     *            NLOOP,ILOOP,NB3,DIVANG,MXANGM,
     *            X(ILBLINT),X(INUMINT))
C
  130 CONTINUE
      CALL RETFMP(NEEDB)
      IF(EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9110)
         CALL TIMIT(1)
      END IF
C
C     ----- ALLOCATE MEMORY FOR -----
C     ----- PIJ, PAB, PAI, WIJ, WAB, WAI, PMN, WMN, AND LAI -----
C
      CALL VALFMP(LOADFM)
      IPIJ = 1    + LOADFM
      IPAB = IPIJ + NOCC2
      IPAI = IPAB + NVIR2
      IWIJ = IPAI + NOCC*NVIR
      IWAB = IWIJ + NOCC2
      IWAI = IWAB + NVIR2
      IPMN = IWAI + NOCC*NVAC
      IWMN = IPMN + NBF2
      ILAI = IWMN + NBF2
      LAST = ILAI + NOCC*NVIR
      NEEDC= LAST - LOADFM -1
      CALL GETFMP(NEEDC)
C
C     --- ALLOCATE MEMORY FOR 4TH TRANSFORMATION AND MP2 AMPLITUDE ---
C     THERE ARE ADDITIONAL ALLOCATIONS INSIDE -MPGQTB-,
C     SO CHECK RUNS CAN'T SIMPLY BRANCH AROUND ITS CALL
C
      CALL VALFMP(LOADFM)
      IAIBJC= 1      + LOADFM
      IAIBJA= IAIBJC + NOCC*NVIR*NACORE
      IKIBJC= IAIBJA + NOCC*NVIR*NOAC
      IKIBJA= IKIBJC + NOCC2*NACORE
      IWIJ1 = IKIBJA + NOCC2*NOAC
      IPIJ1 = IWIJ1  + NOCC*NOAC
      IPMX  = IPIJ1  + NOAC2
      LAST  = IPMX   + NBF*MAX(NOCC,NVIR)
      NEEDD = LAST   - LOADFM -1
      CALL GETFMP(NEEDD)
C
C     ----- OPEN DIRECT ACCESS FILE (WORK9, WORK16, WORK20) -----
C
      CALL MPOPEN(NFT9 ,NOCC2,'MCQD50')
      CALL MPOPEN(NFT16,NOCC2,'MCQD52')
      CALL MPOPEN(NFT20,NOAC2,'MCQD54')
C
C     ----- 4TH TRANSFORMATION, MP2 ENERGY, -----
C     ----- MP2 AMPLITUDE, PIJ, PAB, WIJ1, WAB1, AND WAI1 -----
C
      CALL MOCOPY(X(IVEC),X(IVECO2),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,6)
      CALL MOCOPY(X(IVEC),X(IVECO3),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,1)
      CALL MOCOPY(X(IVEC),X(IVECV2),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,7)
      CALL MOCOPY(X(IVEC),X(IVECV3),NVIR,NBF,NOCC,NVIR,NACORE,
     *            NLOOP(1,3),ILOOP(1,1,2),8)
      CALL MPGQTB(NFT9,NFT15,NFT16,NFT20,X(IVECO2),X(IVECO3),X(IVECV),
     *            X(IVECV2),X(IVECV3),X(IENG),X(IPIJ),X(IWIJ),X(IPAB),
     *            X(IWAB),X(IWAI),X(IPIJ1),X(IWIJ1),X(IPMX),X(IPMN),
     *            X(ILAI),X(IAIBJC),X(IAIBJA),X(IKIBJC),X(IKIBJA),
     *            NBF,NOCC,NVIR,NACORE,NOAC,NVAC,NLOOP,NCACHE,
     *            LENDDI,LRECV,NPROCS,NEEDQTB)
      NEED3 = NEEDA + NEEDC + NEEDD + NEEDQTB
      CALL RETFMP(NEEDD)
C
C     ----- CLOSE DIRECT ACCESS FILE (WORK15) -----
C
      CALL MPCLOS(NFT15,'DELETE')
C
C     ----- PRINT ENERGY DATA -----
C
      IF(GOPARR) CALL DDI_GSUMF(1000,E2,1)
      EHF =ESCF
      EMP2=EHF+E2
      IF(EXETYP.NE.CHECK .OR. EXESAV.EQ.CHECK) THEN
         IF(MASWRK) WRITE(IW,9090) EHF,E2,EMP2
      END IF
C
      IF(EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9120)
         CALL TIMIT(1)
      END IF
C
C     ----- 1ST, 2ND, 3RD QUARTER TRANSFORMATIONS  -----
C
      NA3 = NCACHE/NOCC3 + 1
      NBT2= NCACHE/NOCC2 + 1
      IF(NA3.GT.MAX(NOCC,NVIR)) NA3= MAX(NOCC,NVIR)
      IF(NBT2.GT.NBF) NBT2= NBF
C
      CALL GOTFMP(NPTR2)
      NNPTR2= (NANGM**3)*NBF+NOCC3*NA3+NOAC2*NBF*(1+NANGM)+2*NOCC*NBF
      NPTR2= NPTR2-NNPTR2
      MPTR2= NPTR2-(NANGM**2)*NOCC*NBF
C     IF(MPTR2.LT.0.AND.MASWRK) THEN
      IF(MPTR2.LT.0) THEN
         WRITE(IW,9040) (NEEDA+NEEDC+NNPTR2
     *                 +(NANGM**2)*NOCC*NBF+NOCC3*NBF)/1000000
         CALL ABRT
      ENDIF
      JPTR2= NANGM*NOCC3*NBF
      DIVANG=.FALSE.
      IF(MPTR2.LT.JPTR2) THEN
         DO 210 I=1, NANGM-1
            MXANGM= NANGM-I
            JPTR2= MXANGM*NOCC3*NBF
            IF(MPTR2.GE.JPTR2) GO TO 220
  210    CONTINUE
         IF(MASWRK) WRITE(IW,9040) (NEEDA+NEEDC+JPTR2)/1000000
         CALL ABRT
  220    DIVANG=.TRUE.
         IF(MASWRK) WRITE(IW,9050) MXANGM
      ENDIF
      MPTR2= NPTR2-JPTR2
      NS1= MPTR2/(NANGM*NOCC*NBF)
      IF(NS1.GT.NBF) NS1= NBF
C
      CALL VALFMP(LOADFM)
      IINT  = 1     + LOADFM
      IPTR1 = IINT  +(NANGM**3)*NBF
      IPTR2 = IPTR1 + NANGM*NS1*NOCC*NBF
      IPTR3 = IPTR2 + JPTR2
      ITLB  = IPTR3 + NA3*NOCC3
      ITLS  = ITLB  + NOAC2*NBF
      ILNI1 = ITLS  + NOAC2*NBF*NANGM
      ILNI2 = ILNI1 + NOCC*NBF
      LAST  = ILNI2 + NOCC*NBF
      NEEDE = LAST  - LOADFM -1
      CALL GETFMP(NEEDE)
      NEED4 = NEEDA + NEEDC + NEEDE
C
      IF(EXETYP.EQ.CHECK) GO TO 260
C
C     ----- OPEN DIRECT ACCESS FILE (WORK15, WORK19) -----
C
      CALL MPOPEN(NFT15,NOCC3,'MCQD51')
      CALL MPOPEN(NFT19,NOCC3,'MCQD53')
C
      CALL MOCOPY(X(IVEC),X(IVECO2),NOCC,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,5)
      CALL MOCOPY(X(IVEC),X(IVECV2),NVIR,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,4)
      CALL MPGQTC(NFT15,NFT19,NFT20,X(IVECO),X(IVECO2),X(IVECO3),
     *            X(IVECV),X(IVECV2),X(ILAI),X(ILNI1),X(ILNI2),X(IGHND),
     *            X(IXINT),X(IINT),X(IPTR1),X(IPTR2),X(IPTR3),X(ITLB),
     *            X(ITLS),X(IPMN),NOCC,NVIR,NBF,NOAC,NVAC,NACORE,NANGM,
     *            CUTOFF,NLOOP,ILOOP,NA3,NS1,NBT2,DIVANG,MXANGM,
     *            X(ILBLINT),X(INUMINT))
C
  260 CONTINUE
      CALL RETFMP(NEEDE)
      IF(EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9130)
         CALL TIMIT(1)
      END IF
C
C     --- ALLOCATE MEMORY FOR 4TH TRANSFORMATION (IJ|KL) AND (IJ|AB) ---
C
      NWORK= NCACHE/NOCC2+1
      IF(NWORK.GT.NVIR) NWORK= NVIR
      MAXQ= 1
      IF(GOPARR) THEN
         DO I= 1,NPROC
            MAXQ= MAX(MAXQ,NLOOP(I,4))
         ENDDO
      ENDIF
      CALL VALFMP(LOADFM)
      IPTR3 = 1     + LOADFM
      IDDI  = IPTR3 + NOCC3*NBF
      IPTRO4= IDDI  + NOCC3*MAXQ*NPROC
      IPTRV4= IPTRO4+ NOCC3*NOCC
      IAIBJ = IPTRV4+ NOCC3*NVIR
      IWIJ1 = IAIBJ + NOCC2*NVIR
      IWIJ2 = IWIJ1 + NOCC3
      IWORK = IWIJ2 + NOCC2
      IENGAI= IWORK + NOCC2*NWORK
      LAST  = IENGAI+ NOCC*NVIR
      NEEDF = LAST  - LOADFM -1
      CALL GETFMP(NEEDF)
      NEED5 = NEEDA + NEEDC + NEEDF
C
      IF(EXETYP.EQ.CHECK) GO TO 270
C
C     ----- 4TH TRANSFORMATION (IJ|KL) AND (IJ|AB) -----
C
      CALL MOCOPY(X(IVEC),X(IVECO),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP(1,3),ILOOP(1,1,2),9)
      CALL MOCOPY(X(IVEC),X(IVECV),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP(1,3),ILOOP(1,1,2),10)
      CALL MPGQTD(NFT9,NFT15,NFT19,X(IVECO),X(IVECV),X(IENG),X(IPIJ),
     *            X(IPAB),X(IPTR3),X(IDDI),X(IPTRO4),X(IPTRV4),X(IAIBJ),
     *            X(IWIJ),X(IWIJ1),X(IWIJ2),X(IWORK),X(IENGAI),NOCC,
     *            NVIR,NBF,NOCC3,NLOOP,MAXQ,NWORK,LENDDI,LRECV,NPROCS)
      CALL MPCLOS(NFT15,'DELETE')
      CALL MPCLOS(NFT19,'DELETE')
C
  270 CONTINUE
      CALL RETFMP(NEEDF)
      IF(EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9140)
         CALL TIMIT(1)
      END IF
C
C     ----- ALLOCATE MEMORY FOR SOLVING CPHF EQUATIONS OF MP2 -----
C
      CALL VALFMP(LOADFM)
      IAIBJ  = 1      + LOADFM
      IGMRES1= IAIBJ  + NOCC2*NVIR
      IGMRES2= IGMRES1+ ITERMX*(ITERMX+1)
      IGMVEC = IGMRES2+ ITERMX*(ITERMX+1)
      IGMTRN = IGMVEC + NOCC*NVIR*(ITERMX+2)
      IGMWRK1= IGMTRN +(ITERMX+1)*(ITERMX+1)
      IGMWRK2= IGMWRK1+(ITERMX+1)*(ITERMX+1)
      IGMFAC = IGMWRK2+(ITERMX+1)*(ITERMX+1)
      LAST   = IGMFAC + ITERMX
      NEEDG  = LAST   - LOADFM -1
      CALL GETFMP(NEEDG)
      NEED6 = NEEDA + NEEDC + NEEDG
C
      IF(EXETYP.EQ.CHECK) GO TO 280
C
C     ----- SOLVE CPHF EQUATIONS OF MP2 -----
C
      CALL MPGCPHF(NFT9,X(IENG),X(IPAI),X(ILAI),X(IAIBJ),X(IGMRES1),
     *             X(IGMRES2),X(IGMVEC),X(IGMTRN),X(IGMWRK1),X(IGMWRK2),
     *             X(IGMFAC),NOCC,NVIR,ITERMX)
C
  280 CONTINUE
      CALL RETFMP(NEEDG)
C
C     ----- ALLOCATE MEMORY FOR WIJ[III] AND WAI[II] -----
C
      CALL VALFMP(LOADFM)
      IWIJ1 = 1     + LOADFM
      IXKIBJ= IWIJ1 + NOCC2
      LAST  = IXKIBJ+ NOCC2*NOCC
      NEEDH = LAST  - LOADFM -1
      CALL GETFMP(NEEDH)
      NEED7 = NEEDA + NEEDC + NEEDH
C
      IF(EXETYP.EQ.CHECK) GO TO 290
C
C     ----- CALCULATE WIJ[III] AND WAI[II] -----
C
      CALL MPGWAI(NFT16,X(IENG),X(IPAI),X(IWIJ),X(IWAI),X(IWIJ1),
     *            X(IXKIBJ),NOCC,NVIR)
      CALL MPCLOS(NFT16,'DELETE')
      CALL MPCLOS(NFT9,'DELETE')
C
  290 CONTINUE
      CALL RETFMP(NEEDH)
C
C     ----- ALLOCATE MEMORY FOR WMN AND PMN -----
C
      CALL VALFMP(LOADFM)
      IWMJ = 1    + LOADFM
      IWMB = IWMJ + NOCC*NBF
      IPMJ = IWMB + NVIR*NBF
      IPMB = IPMJ + NOCC*NBF
      IPAO = IPMB + NVIR*NBF
      IWAO = IPAO + NBF3
      LAST = IWAO + NBF3
      NEEDI= LAST - LOADFM -1
      CALL GETFMP(NEEDI)
      NEED8 = NEEDA + NEEDC + NEEDI
C
      IF(EXETYP.EQ.CHECK) GO TO 300
C
C     ----- CALCULATE WMN AND PMN -----
C
      CALL MOCOPY(X(IVEC),X(IVECO),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,1)
      CALL MOCOPY(X(IVEC),X(IVECO2),NOCC,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,5)
      CALL MOCOPY(X(IVEC),X(IVECV),NBF,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,2)
      CALL MPGTAO(X(IVECO),X(IVECO2),X(IVECV),X(IENG),X(IPMN),X(IWMN),
     *             X(IPAI),X(IWIJ),X(IWAB),X(IWAI),X(IWMJ),X(IWMB),
     *             X(IPAO),X(IWAO),NOCC,NVIR,NBF)
C
  300 CONTINUE
      CALL RETFMP(NEEDI)
      IF(EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9150)
         CALL TIMIT(1)
      END IF
C
C     ------ CALCULATE GRADIENT OF ONE ELECTRON TERMS -----
C
      IF(ICHKPASS.EQ.2  .OR.  EXESAV.EQ.CHECK) CALL STVDER
C
C     ------ ALLOCATE MEMORY FOR GRADIENT OF TWO ELECTRON TERMS -----
C     THERE ARE ADDITIONAL ALLOCATIONS INSIDE -MPGJKD-,
C     SO CHECK RUNS CAN'T SIMPLY BRANCH AROUND ITS CALL
C
      CALL VALFMP(LOADFM)
      ITIJ = 1    + LOADFM
      ITMJ = ITIJ + NOAC2*NANGM*NANGM
      ITMN = ITMJ + NOAC*NBF*NANGM*NANGM
      IPHF = ITMN + NANGM**4
      IPMP2= IPHF + NBF3
      IDE1 = IPMP2+ NBF3
      LAST = IDE1 + NAT*3
      NEEDJ= LAST - LOADFM -1
      CALL GETFMP(NEEDJ)
C
C     ------ CALCULATE GRADIENT OF TWO ELECTRON TERMS -----
C
      CALL MOCOPY(X(IVEC),X(IVECO),NOCC,NBF,NOCC,NVIR,NACORE,
     *            NLOOP,ILOOP,3)
      CALL MPGJKD(NFT20,X(IVECO),X(IXINT),X(ITIJ),X(ITMJ),X(ITMN),
     *            X(IPHF),X(IPMP2),X(IDE1),NLOOP,ILOOP,NOCC,NBF,NOAC,
     *            NANGM,NEEDGRD)
      CALL MPCLOS(NFT20,'DELETE')
      NEED9 = NEEDA + NEEDC + NEEDJ + NEEDGRD
C
      CALL RETFMP(NEEDJ)
      CALL RETFMP(NEEDC)
      CALL RETFMP(NEEDA)
C        RYS QUADRATURE USES AN EXTRA 49*30**2 ARRAY NOT INCLUDED HERE.
C        THE TINY MEMORY ALLOCATED BY MP2GRDIMS IS NOT INCLUDED HERE.
      MXNEED = MAX(NEED2,NEED3,NEED4,NEED5,NEED6,NEED7,NEED8,NEED9)
      IF(MASWRK) WRITE(IW,9025) MXNEED,(MXNEED-1)/1000000+1
C
C        REAL RUNS NOW BRANCH BACK TO DO THE ACTUAL CALCULATION
C
      EXETYP=EXESAV
      IF(ICHKPASS.EQ.1  .AND.  EXESAV.NE.CHECK) GO TO 50
C
      DSKWRK=DSKSAV
      CALL SYMON
C
      RETURN
C
 9000 FORMAT(/1X,31('-'),3X,29(1H-)/
     *        1X,'MP2 ENERGY GRADIENT CALCULATION',3X,
     *           'PROGRAM WRITTEN BY K.ISHIMURA'/
     *        1X,31('-'),3X,29(1H-)/
     *        1X,'NUMBER OF OCCUPIED ORBITALS            =',I5,/
     *        1X,'NUMBER OF VIRTUAL ORBITALS             =',I5,/
     *        1X,'NUMBER OF ACTIVE OCCUPIED ORBITALS     =',I5,/
     *        1X,'NUMBER OF ACTIVE VIRTUAL ORBITALS      =',I5,/
     *        1X,'NUMBER OF BASIS FUNCTIONS              =',I5,/
     *        1X,'NUMBER OF BASIS SHELLS                 =',I5,/
     *        1X,'INTEGRAL THRESHOLD                     =',E8.1)
C9005 FORMAT(/1X,'NUMBER OF REDUNDANT MOS IN FMO         =',I5)
C9010 FORMAT(/1X,'MEMORY AND DISK REQUIREMENTS FOR SEMI-DIRECT MP2',1X,
C    *        'TRANSFORMATION'/
C    *        '  MAXIMUM MEMORY=',I10,' WORDS (',I7,' MWORDS)'/
C    *        '  MINIMUM MEMORY=',I10,' WORDS (',I7,' MWORDS)')
C9020 FORMAT(/1X,'AVAILABLE MEMORY =', I10,' MWORDS'/
C    *        1X,'NOT ENOUGH MEMORY FOR MP2 GRADIENT'/)
 9015 FORMAT(1X,'  REQUIRED DISK =',I10,' MBYTES/PROCESSOR'/)
 9025 FORMAT(1X,'MAXIMUM MEMORY USED DURING MP2 GRADIENT=',I12,
     *           ' WORDS, OR',I8,' MWORDS.')
 9030 FORMAT(/1X,'*** WARNING: THE MAXIMUM LCAO COEF=',F15.6,
     *           ', FOR AO',I6,' IN MO=',I6,','/
     *        1X,'EXCEEDS 100.  ADDITIONAL INPUTS MAY BE REQUIRED TO',
     *           ' OBTAIN GOOD MP2 ACCURACY:'/
     *        1X,' $CONTRL ICUT=10 OR 11'/
     *        1X,' $SCF    CONV=2.5D-7 FDIFF=.FALSE.'/)
 9040 FORMAT(1X,'**** ERROR: TOO LITTLE MEMORY TO RUN MP2 ',
     *          'GRADIENT *****'/
     *       1X,'INCREASE MEMORY TO AT LEAST',I7,'MWORDS.')
 9050 FORMAT(1X,'MULTIPLE PASSES',
     *      /1X,'MAXIMUM ANGULAR MOMENTUM=',I3)
 9090 FORMAT(/1X,'RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE'/
     *       13X,'  E(0)=',1X,F20.10/
     *       13X,'  E(1)=',9X,'0.0'/
     *       13X,'  E(2)=',1X,F20.10/
     *       13X,'E(MP2)=',1X,F20.10)
 9110 FORMAT(1X,'.... DONE WITH 1ST THREE INDEX TRANSFORMATIONS ....')
 9120 FORMAT(1X,'.... DONE WITH 4TH INDEX TRANSFORMATION AND',
     *          ' ENERGY ....')
 9130 FORMAT(1X,'.... DONE WITH 1ST THREE INDEX TRANSFORMATIONS FOR',
     *          ' GRADIENT ....')
 9140 FORMAT(1X,'.... DONE WITH 4TH INDEX TRANSFORMATION FOR',
     *          ' GRADIENT ....')
 9150 FORMAT(1X,'.... DONE WITH RESPONSE EQUATIONS AND DENSITY ....')
C
      END
C
C*MODULE MP2GR2  *DECK MOCOPY
      SUBROUTINE MOCOPY(VEC,VECX,IDIM,NBF,NOCC,NVIR,NACORE,
     *                  NLOOP,ILOOP,ITYPE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,PARDYN
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      DIMENSION VEC(NBF,*),VECX(IDIM,*),NLOOP(NPROC),
     *          ILOOP(NSHELL,NPROC)
C
      PARDYN= GOPARR.AND.IBTYP.EQ.1
C
C     COPY OCCUPIED MO COEFFICIENTS
C
      IF((ITYPE.EQ.1).OR.(ITYPE.EQ.6.AND..NOT.PARDYN).OR.
     *   (ITYPE.EQ.9.AND..NOT.GOPARR)) THEN
        CALL DCOPY(NOCC*NBF,VEC,1,VECX,1)
C
C     COYP VIRTUAL MO COEFFICIENTS
C
      ELSEIF((ITYPE.EQ.2).OR.(ITYPE.EQ.7.AND..NOT.PARDYN).OR.
     *       (ITYPE.EQ.10.AND..NOT.GOPARR)) THEN
        CALL DCOPY(NVIR*NBF,VEC(1,NOCC+1),1,VECX,1)
C
C     TRANSPOSE OCCUPIED MO COEFFICIENTS (ORDER:ACTIVE,CORE)
C
      ELSEIF(ITYPE.EQ.3) THEN
        NOAC= NOCC-NACORE
        DO I= 1,NACORE
          DO J= 1,NBF
            VECX(NOAC+I,J)= VEC(J,I)
          ENDDO
        ENDDO
        DO I= 1,NOAC
          DO J= 1,NBF
            VECX(I,J)= VEC(J,I+NACORE)
          ENDDO
        ENDDO
C
C     TRANSPOSE VIRTUAL MO COEFFICIENTS
C
      ELSEIF((ITYPE.EQ.4).OR.(ITYPE.EQ.8.AND..NOT.GOPARR)) THEN
        DO I= 1,NVIR
          DO J= 1,NBF
            VECX(I,J)= VEC(J,I+NOCC)
          ENDDO
        ENDDO
C
C     TRANSPOSE OCCUPIED MO COEFFICIENTS
C
      ELSEIF(ITYPE.EQ.5) THEN
        DO I= 1,NOCC
          DO J= 1,NBF
            VECX(I,J)= VEC(J,I)
          ENDDO
        ENDDO
C
C     SORT OCCUPIED MO COEFFICIENTS FOR 4TH TRANSFORMATION
C
      ELSEIF((ITYPE.EQ.6.AND.PARDYN).OR.
     *       (ITYPE.EQ.9.AND.GOPARR)) THEN
        N= 0
        DO I= 1,NPROC
          DO J= 1,NLOOP(I)
            ISH= ILOOP(J,I)
            MINI= KMIN(ISH)
            MAXI= KMAX(ISH)
            LOCI= KLOC(ISH)-MINI
            DO K= MINI,MAXI
              N= N+1
              DO L= 1,NOCC
                VECX(N,L)= VEC(K+LOCI,L)
              ENDDO
            ENDDO
          ENDDO
        ENDDO
C
C     SORT VIRTUAL MO COEFFICIENTS FOR 4TH TRANSFORMATION
C
      ELSEIF((ITYPE.EQ.7.AND.PARDYN).OR.
     *       (ITYPE.EQ.10.AND.GOPARR)) THEN
        N= 0
        DO I= 1,NPROC
          DO J= 1,NLOOP(I)
            ISH= ILOOP(J,I)
            MINI= KMIN(ISH)
            MAXI= KMAX(ISH)
            LOCI= KLOC(ISH)-MINI
            DO K= MINI,MAXI
              N= N+1
              DO L= 1,NVIR
                VECX(N,L)= VEC(K+LOCI,L+NOCC)
              ENDDO
            ENDDO
          ENDDO
        ENDDO
C
C     SORT VIRTUAL MO COEEFICIENTS FOR BACK-TRANSFORMATION
C
      ELSEIF(ITYPE.EQ.8.AND.GOPARR) THEN
         N= 0
         DO I= 1,NPROC
            DO J= 1,NLOOP(I)
               ISH= ILOOP(J,I)
               MINI= KMIN(ISH)
               MAXI= KMAX(ISH)
               LOCI= KLOC(ISH)-MINI
               DO K= MINI,MAXI
                  N= N+1
                  DO L= 1,NVIR
                     VECX(L,N)= VEC(K+LOCI,L+NOCC)
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGQTA
      SUBROUTINE MPGQTA(NFT15,VECO,VECV,GHONDO,XINTS,
     *                  AOINTS,PTR1,PTR2,PTR3,
     *                  NOCC,NVIR,NBF,NOAC,NANGM,CUTOFF,NUMP,
     *                  NLOOP,ILOOP,NB3,DIVANG,MXANGM,
     *                  LBLINT,NUMINT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SCHSKP
      LOGICAL DIVANG,PASSCH,NXT
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (ZER=0.0D+00)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SHLEXC/ NORGSH(3),NORGSP(3),IEXCH,MANGM,NGTH(4)
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
C
      DIMENSION VECO(NOCC,*),VECV(NBF,*),GHONDO(*),XINTS(*),
     *          AOINTS(*),PTR1(*),PTR2(*),PTR3(*),NUMP(0:NPROC)
      DIMENSION NLOOP(0:NPROC-1,4),ILOOP(NSHELL,0:NPROC-1,2)
      DIMENSION LBLINT(NBF*(NANGM**3)),NUMINT(NANGM**3)
C
C     1ST, 2ND, 3RD QUARTER TRANSFORMATIONS
C
      NXT = IBTYP.EQ.1
      NEXT = -1
      MINE = -1
C
      NORG  = 0
      IEXCH = 1
      ISTART= 1
      IEND  = NSHELL
      IREC  = 0
C
      TAO=ZER
      TQT1=ZER
      TQT2=ZER
      TQT3=ZER
C
C     ----- STATIC PARALLEL -----
C
      IF((.NOT.NXT) .AND. GOPARR) THEN
         ISTART= NUMP(ME)+1
         IEND  = NUMP(ME+1)
      ENDIF
C
C     ----- DYNAMIC PARALLEL -----
C
      IF(NXT .AND. GOPARR) THEN
         NLOOP(ME,1)= 0
         NLOOP(ME,2)= 0
      ENDIF
C
C
C     ----- START I SHELL -----
C
      DO 500 ISH= ISTART,IEND
C
         IF(MOD(ISH,25).EQ.1) THEN
            IF(MASWRK) WRITE(IW,9999) ISH
            CALL FLSHBF(IW)
         END IF
C
C     ----- GO PARALLEL! -----
C
         IF(NXT.AND.GOPARR) THEN
            MINE = MINE + 1
            IF (MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF (NEXT.NE.MINE) GO TO 500
            NLOOP(ME,1)= NLOOP(ME,1)+1
            ILOOP(NLOOP(ME,1),ME,1)=ISH
         END IF
C
         MINI= KMIN(ISH)
         MAXI= KMAX(ISH)
         MINI2= MINI
         MAXI2= MAXI
         ITMP= MAXI2-MINI2+1
         IDIV= 0
  110    IF(DIVANG.AND.ITMP.GT.MXANGM) THEN
            MINI2= MXANGM*IDIV+KMIN(ISH)
            MAXI2= MXANGM+MINI2-1
            IF(MAXI2.GT.KMAX(ISH)) MAXI2=KMAX(ISH)
            IDIV= IDIV+1
         ENDIF
         NUMI= MAXI2-MINI2+1
         IF(NXT.AND.GOPARR) NLOOP(ME,2)=NLOOP(ME,2)+NUMI
         CALL VCLR(PTR2,1,NOCC*NOCC*NBF*NUMI)
C
C     ----- START K SHELL -----
C
      DO 400 KSH= 1,NSHELL
         MINK= KMIN(KSH)
         MAXK= KMAX(KSH)
         NUMK= MAXK-MINK+1
         LOCK= KLOC(KSH)-MINK
C
C     ----- START L SHELL -----
C
      DO 300 LSH= 1,KSH
         MINL= KMIN(LSH)
         MAXL= KMAX(LSH)
         NUML= MAXL-MINL+1
         LOCL= KLOC(LSH)-MINL
         KLKL= IA(KSH)+LSH
         PASSCH=.FALSE.
         DO II=1,NUMI*NUMK*NUML
            NUMINT(II)=0
         ENDDO
C
C     ----- START J SHELL -----
C
      DO 200 JSH= 1,NSHELL
         IJIJ= IA(ISH)+JSH
         IF(JSH.GT.ISH) IJIJ= IA(JSH)+ISH
         TEST= XINTS(IJIJ)*XINTS(KLKL)
         SCHSKP= TEST.LT.CUTOFF
         MINJ= KMIN(JSH)
         MAXJ= KMAX(JSH)
         LOCJ= KLOC(JSH)-MINJ
C
         IF(.NOT.SCHSKP) THEN
            PASSCH= .TRUE.
CTIME
            CALL TSECND(TAO0)
            CALL SHELLQUARTMP2(GHONDO)
C
C     ----- SAVE AO INTEGRALS -----
C
            CALL SAVEAO(AOINTS,GHONDO,CUTOFF,
     *                  LBLINT,NUMINT,KSH,LSH,LOCJ,MINI,
     *                  MINI2,MAXI2,MINJ,MAXJ,MINK,MAXK,MINL,MAXL)
CTIME
            CALL TSECND(TAO1)
            TAO=TAO+(TAO1-TAO0)
C
         ENDIF
  200 CONTINUE
C
C     ----- END J SHELL -----
C
      IF(.NOT.PASSCH) GO TO 300
C
C     ----- 1ST QUARTER TRANSFORMATION -----
C
CTIME
      CALL TSECND(TQT1S)
      CALL MPGT1A(VECO,AOINTS,PTR1,NOCC,NBF,NUMI,NUMK,NUML,
     *            NANGM,LBLINT,NUMINT,KSH,LSH)
CTIME
      CALL TSECND(TQT1E)
      TQT1=TQT1+(TQT1E-TQT1S)
C
C     ----- 2ND QUARTER TRANSFORMATION -----
C
CTIME
      CALL TSECND(TQT2S)
      CALL MPGT2A(VECO,PTR1,PTR2,
     *            NOCC,NBF,NUMI,NUMK,NUML,NANGM,CUTOFF,
     *            KSH,LSH,MINK,MINL,LOCK,LOCL)
CTIME
      CALL TSECND(TQT2E)
      TQT2=TQT2+(TQT2E-TQT2S)
C
  300 CONTINUE
C
C     ----- END L SHELL -----
C
  400 CONTINUE
C
C     ----- END K SHELL -----
C
C
C     ----- 3RD QUARTER TRANSFORMATION -----
C
CTIME
      CALL TSECND(TQT3S)
      CALL MPGT3A(NFT15,VECV,PTR2,PTR3,NOCC,NVIR,NBF,NOAC,IREC,NUMI,
     *            NB3)
CTIME
      CALL TSECND(TQT3E)
      TQT3=TQT3+(TQT3E-TQT3S)
C
      IF(DIVANG.AND.MAXI2.NE.KMAX(ISH)) GO TO 110
C
  500 CONTINUE
C
C     ----- END I SHELL -----
C
C     ----- CHECK DYNAMIC PARALLEL INFORMATIONS -----
C
      IF(NXT.AND.GOPARR) THEN
         CALL DDI_DLBRESET
         ISIZE= 8/NWDVAR
         ILSIZE= NSHELL*ISIZE
         DO 600 I=1,NPROC-1
            IF(MASWRK) THEN
               CALL DDI_RECV(NLOOP(I,1),ISIZE,I)
               CALL DDI_RECV(NLOOP(I,2),ISIZE,I)
               CALL DDI_RECV(ILOOP(1,I,1),ILSIZE,I)
            ELSE
               IF(I.EQ.ME)CALL DDI_SEND(NLOOP(I,1),ISIZE,MASTER)
               IF(I.EQ.ME)CALL DDI_SEND(NLOOP(I,2),ISIZE,MASTER)
               IF(I.EQ.ME)CALL DDI_SEND(ILOOP(1,I,1),ILSIZE,MASTER)
            ENDIF
  600    CONTINUE
         CALL DDI_BCAST(1012,'I',NLOOP,4*NPROC,MASTER)
         CALL DDI_BCAST(1011,'I',ILOOP,NSHELL*NPROC,MASTER)
      ENDIF
C
CTIME
CTIME WRITE(IW,999)TAO,TQT1,TQT2,TQT3
C 999 FORMAT(/1X,'TIME TO GENERATE AO INTEGRALS=     ',F9.1,
C    *       /1X,'TIME OF 1ST QUARTER TRANSFORMATION=',F9.1,
C    *       /1X,'TIME OF 2ND QUARTER TRANSFORMATION=',F9.1,
C    *       /1X,'TIME OF 3RD QUARTER TRANSFORMATION=',F9.1)
 9999 FORMAT(1X,I5,'ST SHELL LOOP START')
      RETURN
      END
C
C
C*MODULE MP2GR2  *DECK MPGT1A
      SUBROUTINE MPGT1A(VECO,AOINTS,PTR1,NOCC,NBF,NUMI,NUMK,NUML,
     *                  NANGM,LBLINT,NUMINT,KSH,LSH)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZER=0.0D+00)
C
      LOGICAL KANDL
C
      DIMENSION VECO(NOCC,*),AOINTS(NBF,*),PTR1(NOCC,NANGM,NANGM,NANGM)
      DIMENSION LBLINT(NBF,*),NUMINT(*)
C
C     1ST QUARTER TRANSFORMATION
C
      KANDL= KSH.EQ.LSH
      LNUM = NUML
C
      IKL= 0
      DO 160 I= 1,NUMI
         DO 150 K= 1,NUMK
            IF(KANDL) LNUM= K
            DO 140 L= 1,LNUM
               IKL= IKL+1
               DO 110 MOI= 1,NOCC
                  PTR1(MOI,L,K,I)= ZER
  110          CONTINUE
               DO 130 J= 1,NUMINT(IKL)
                  JLABEL= LBLINT(J,IKL)
                  DUM   = AOINTS(J,IKL)
                  DO 120 MOI= 1,NOCC
                     PTR1(MOI,L,K,I)= PTR1(MOI,L,K,I)
     *                               +DUM*VECO(MOI,JLABEL)
  120             CONTINUE
  130          CONTINUE
  140       CONTINUE
  150    CONTINUE
  160 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT2A
      SUBROUTINE MPGT2A(VECO,PTR1,PTR2,
     *                  NOCC,NBF,NUMI,NUMK,NUML,NANGM,CUTOFF,
     *                  KSH,LSH,MINK,MINL,LOCK,LOCL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL KNOTL
C
      DIMENSION VECO(NOCC,*),PTR1(NOCC,NANGM,NANGM,NANGM),
     *          PTR2(NOCC,NOCC,NBF,*)
C
C     2ND QUARTER TRANSFORMATION
C
      KNOTL= KSH.NE.LSH
      LOCK2= LOCK+MINK-1
      LOCL2= LOCL+MINL-1
C
C
      IF(KNOTL)THEN
      DO 150 I= 1,NUMI
         DO 140 L= 1,NUML
            LL= LOCL2+L
            DO 130 K= 1,NUMK
               KK= LOCK2+K
               DO 120 MOI= 1,NOCC
                  DUM1= PTR1(MOI,L,K,I)
CISHI             IF(ABS(DUM1)*TMP1.LT.CUTOFF) GO TO 120
                  IF(ABS(DUM1).LT.CUTOFF) GO TO 120
                  DO 110 MOJ= 1,NOCC
                     PTR2(MOJ,MOI,LL,I)= PTR2(MOJ,MOI,LL,I)
     *                                  +DUM1*VECO(MOJ,KK)
  110             CONTINUE
  120          CONTINUE
  130       CONTINUE
  140    CONTINUE
  150 CONTINUE
C
      DO 250 I=1,NUMI
         DO 240 K=1,NUMK
            KK= LOCK2+K
            DO 230 L=1,NUML
               LL= LOCL2+L
               DO 220 MOI=1,NOCC
                  DUM1= PTR1(MOI,L,K,I)
CISHI             IF(ABS(DUM1)*TMP1.LT.CUTOFF) GO TO 220
                  IF(ABS(DUM1).LT.CUTOFF) GO TO 220
                  DO 210 MOJ= 1,NOCC
                     PTR2(MOJ,MOI,KK,I)= PTR2(MOJ,MOI,KK,I)
     *                                  +DUM1*VECO(MOJ,LL)
  210             CONTINUE
  220          CONTINUE
  230       CONTINUE
  240    CONTINUE
  250 CONTINUE
C
      ELSE
      DO 370 I=1,NUMI
         DO 360 K= 1,NUMK
            LMAX= K
            KK= LOCK2+K
            DO 350 L= 1,LMAX
               LL= LOCL2+L
               DO 340 MOI= 1,NOCC
                  DUM1= PTR1(MOI,L,K,I)
                  DUM3= ABS(DUM1)
CISHI             IF(DUM3*TMP1.LT.CUTOFF) GO TO 320
                  IF(DUM3.LT.CUTOFF) GO TO 320
                  DO 310 MOJ= 1,NOCC
                     PTR2(MOJ,MOI,LL,I)= PTR2(MOJ,MOI,LL,I)
     *                                  +DUM1*VECO(MOJ,KK)
  310             CONTINUE
CI320             IF(DUM3*TMP2.LT.CUTOFF) GO TO 340
  320             IF(DUM3.LT.CUTOFF) GO TO 340
                  DO 330 MOJ= 1,NOCC
                     PTR2(MOJ,MOI,KK,I)= PTR2(MOJ,MOI,KK,I)
     *                                  +DUM1*VECO(MOJ,LL)
  330             CONTINUE
  340          CONTINUE
  350       CONTINUE
  360    CONTINUE
  370 CONTINUE
      ENDIF
C
      RETURN
      END
C
C
C*MODULE MP2GR2  *DECK MPGT3A
      SUBROUTINE MPGT3A(NFT15,VECV,PTR2,PTR3,NOCC,NVIR,NBF,NOAC,IREC,
     *                  NUMI,NB3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER(ZER=0.0D+00, ONE=1.0D+00)
C
      DIMENSION VECV(NBF,*),PTR2(NOCC*NOCC*NBF,*),PTR3(NOCC,NOCC,NB3)
C
C     3RD QUARTER TRANSFORMATION
C
      NOCC2=NOCC*NOCC
      NPASS= NVIR/NB3
      MODB = MOD(NVIR,NB3)
      IF(MODB.NE.0) THEN
         NPASS= NPASS+1
      ELSE
         MODB= NB3
      ENDIF
C
      DO 130 I= 1,NUMI
         MB3= NB3
         DO 120 IPASS= 1,NPASS
            IF(IPASS.EQ.NPASS) MB3= MODB
            MOB= (IPASS-1)*NB3+1
            CALL DGEMM('N','N',NOCC2,MB3,NBF,ONE,PTR2(1,I),NOCC2,
     *                 VECV(1,MOB),NBF,ZER,PTR3,NOCC2)
            DO 110 IWRITE= 1,MB3
               IREC= IREC+1
               WRITE(UNIT=NFT15,REC=IREC)
     *              ((PTR3(MOJ,MOI,IWRITE),MOJ=1,NOAC),MOI=1,NOCC),
     *              ((PTR3(MOJ,MOI,IWRITE),MOJ=NOAC+1,NOCC),MOI=1,NOCC)
  110       CONTINUE
  120    CONTINUE
  130 CONTINUE
C
      RETURN
      END
C
C
C*MODULE MP2GR2  *DECK MPGQTB
      SUBROUTINE MPGQTB(NFT9,NFT15,NFT16,NFT20,VECO,VECO2,VECV,VECV2,
     *                  VECV3,ENG,PIJ,WIJ,PAB,WAB,WAI,
     *                  PIJ1,WIJ1,PMX,PMN,XLAI,
     *                  XAIBJC,XAIBJA,XKIBJC,XKIBJA,
     *                  NBF,NOCC,NVIR,NACORE,NOAC,NVAC,NLOOP,NCACHE,
     *                  LENDDI,LRECV,NPROCS,NEEDQTB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
C     PARAMETER (MXRT=100)
      PARAMETER (ONE=1.0D+00, HALF=0.5D+00,ZER=0.0D+00)
C
C     COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
C    *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FMCOM / X(1)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DIMENSION VECO(NBF,*),VECO2(NBF,*),VECV(NBF,*),VECV2(NBF,*),
     *          VECV3(NVIR,*),ENG(*),PIJ(NOCC,*),WIJ(NOCC,*),
     *          PAB(NVIR,*),WAB(NVIR,*),WAI(NOCC,*),PIJ1(NOAC,*),
     *          WIJ1(NOAC,*),PMX(NBF,*),PMN(NBF,*),XLAI(*),XAIBJC(*),
     *          XAIBJA(*),XKIBJC(*),XKIBJA(*),NLOOP(0:NPROC-1,4)
      DIMENSION LENDDI(0:NPROCS-1,3),LRECV(0:NPROCS-1,2)
C
      DATA CHECK/8HCHECK   /
C
C     4TH QUARTER TRANSFORMATION AND MP2 ENERGY CALCULATION
C
      IRECO=0
      IRECV=0
      MAXP=1
      MAXQ=1
C
      CALL VCLR(PIJ,1,NOCC*NOCC)
      CALL VCLR(WIJ,1,NOCC*NOCC)
      CALL VCLR(PAB,1,NVIR*NVIR)
      CALL VCLR(WAB,1,NVIR*NVIR)
      CALL VCLR(WAI,1,NOCC*NVAC)
      CALL VCLR(PIJ1,1,NOAC*NOAC)
      CALL VCLR(WIJ1,1,NOAC*NOCC)
      CALL VCLR(XLAI,1,NOCC*NVIR)
C
C     ----- CALCULATE SIZE AND POINTER FOR DDI_ISEND AND DDI_IRECV -----
C
      IF(GOPARR) THEN
         DO I= 0,NPROC-1
            MAXP= MAX(MAXP,NLOOP(I,2))
            MAXQ= MAX(MAXQ,NLOOP(I,4))
         ENDDO
         NUMI = NLOOP(ME,2)
         NUMS = NLOOP(ME,4)
         NSUM1= 1
         NSUM2= 1
         DO IPROC= 0,NPROC-1
            LENDDI(IPROC,1)= NLOOP(IPROC,2)*NOCC*NACORE
            LENDDI(IPROC,2)= NLOOP(IPROC,2)*NOCC*NOAC
            LENDDI(IPROC,3)= NLOOP(IPROC,4)*NOAC*NOAC
            LRECV(IPROC,1) = NSUM1
            LRECV(IPROC,2) = NSUM2
            NSUM1          = NSUM1+NLOOP(IPROC,2)
            NSUM2          = NSUM2+NLOOP(IPROC,4)
         ENDDO
      ENDIF
C
      NUMX= NCACHE/(NOCC*NOCC)+1
C
      NUMV= NUMX
      IF(NUMX.GT.NVIR) NUMV= NVIR
      NCYCV= NVIR/NUMV
      MODV = MOD(NVIR,NUMV)
      IF(MODV.NE.0) THEN
         NCYCV= NCYCV+1
      ELSE
         MODV = NUMV
      ENDIF
C
      NUMO= NUMX
      IF(NUMX.GT.NOCC) NUMO= NOCC
      NCYCO= NOCC/NUMO
      MODO = MOD(NOCC,NUMO)
      IF(MODO.NE.0) THEN
         NCYCO= NCYCO+1
      ELSE
         MODO = NUMO
      ENDIF
C
      NPASS = NVIR/NPROC
      IF(MOD(NVIR,NPROC).NE.0) NPASS= NPASS+1
C
      DO IPASS= 1,NPASS
C
C     ----- ALLOCATE MEMORY FOR 4TH TRANSFORMATION -----
C
         CALL VALFMP(LOADFM)
         IDDIC  = 1      + LOADFM
         IDDIA  = IDDIC  + NOCC*NPROC*MAXP*NACORE
         IPTR3C = IDDIA  + NOCC*NPROC*MAXP*NOAC
         IPTR3A = IPTR3C + NOCC*NBF*NACORE
         LAST   = IPTR3A + NOCC*NBF*NOAC
         NEEDA  = LAST   - LOADFM - 1
C             CHECK RUNS CAN ESTIMATE THE MEMORY ABOVE BY NOTING THAT
C             NPROC*MAXP IS APPROXIMATELY NBF, AND NACORE+NOAC = NOCC.
         IF(EXETYP.EQ.CHECK) NEEDA= 2*NOCC*NBF*NOCC
         CALL GETFMP(NEEDA)
         IF(EXETYP.EQ.CHECK) GO TO 210
C
C     ----- 4TH QUARTER TRANSFORMATION -----
C
         MPROC= NPROC
         IF(IPASS.EQ.NPASS) THEN
           MPROC= MOD(NVIR,NPROC)
           IF(MPROC.EQ.0) MPROC= NPROC
         ENDIF
C
         CALL MPGT4A(NFT9,NFT15,NFT16,VECO,VECV2,XAIBJC,XAIBJA,XKIBJC,
     *               XKIBJA,X(IDDIC),X(IDDIA),X(IPTR3C),X(IPTR3A),
     *               LENDDI,LRECV,NOCC,NVIR,NBF,NACORE,NOAC,MAXP,IPASS,
     *               MPROC,NUMI,NUMV,NCYCV,MODV,NUMO,NCYCO,MODO,IRECO,
     *               IRECV)
C
  210    CONTINUE
         CALL RETFMP(NEEDA)
C
C     ----- ALLOCATE MEMORY FOR MP2 ENERGY AND AMPLITUDE ------
C
         CALL VALFMP(LOADFM)
         ITAB  = 1     + LOADFM
         ITLB  = ITAB  + NOAC*NOAC*NVAC
         ITRECV= ITLB  + NOAC*NOAC*NBF
         IDAIBJ= ITRECV+ NOAC*NOAC*MAXQ*NPROC
         LAST  = IDAIBJ+ NOAC*NOAC*NVAC
         NEEDB = LAST  - LOADFM -1
         CALL GETFMP(NEEDB)
         IF(EXETYP.EQ.CHECK) GO TO 220
C
         CALL MPGGEN(ENG,XAIBJA,XKIBJA,PAB,WAB,WAI,PIJ1,WIJ1,
     *               X(ITAB),X(IDAIBJ),
     *               NOCC,NVIR,NACORE,NOAC,NVAC,IPASS)
C
         CALL MPGBT1(NFT20,VECV3,X(ITAB),X(ITLB),X(ITRECV),LENDDI,LRECV,
     *               NVIR,NBF,NOAC,NVAC,MAXQ,IPASS,NUMS,MPROC)
C
  220    CONTINUE
         CALL RETFMP(NEEDB)
         IF(EXETYP.EQ.CHECK) GO TO 300
C
      ENDDO
C
  300 CONTINUE
      NEEDQTB = MAX(NEEDA,NEEDB)
      IF(EXETYP.EQ.CHECK) RETURN
C
C     ----- COPY PIJ, WIJ[I], PAB, AND WAB[I] -----
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(1001,PIJ1,NOAC*NOAC)
         CALL DDI_GSUMF(1002,WIJ1,NOAC*NOCC)
         CALL DDI_GSUMF(1003,PAB,NVIR*NVIR)
         CALL DDI_GSUMF(1004,WAB,NVIR*NVIR)
         CALL DDI_GSUMF(1005,WAI,NOCC*NVAC)
      ENDIF
C
      DO I= 1,NOAC
         DO J= 1,NOAC
            PIJ(J+NACORE,I+NACORE)= PIJ1(J,I)
         ENDDO
      ENDDO
      DO I= 1,NACORE
         DO J= 1,NOAC
            EIJ =-HALF/(ENG(I)-ENG(J+NACORE))
            WIJ(J+NACORE,I)= WIJ1(J,I+NOAC)
            PIJ(J+NACORE,I)= WIJ1(J,I+NOAC)*EIJ
            PIJ(I,J+NACORE)= WIJ1(J,I+NOAC)*EIJ
         ENDDO
      ENDDO
      DO I= 1,NOAC
         DO J= 1,NOAC
            WIJ(J+NACORE,I+NACORE)= WIJ1(J,I)
         ENDDO
      ENDDO
      DO I= NVAC+1,NVIR
         DO J= 1,NVAC
            EAB=-HALF/(ENG(I+NOCC)+ENG(J+NOCC))
            PAB(J,I)= WAB(J,I)*EAB
            PAB(I,J)= WAB(J,I)*EAB
         ENDDO
      ENDDO
C
C
C     ----- CALCULATE WIJ[II] AND WAB[II] -----
C
      DO I= 1,NOCC
         DO J= 1,NOCC
            WIJ(J,I)= WIJ(J,I)+PIJ(J,I)*(ENG(I)+ENG(J))
         ENDDO
      ENDDO
      DO I= 1,NVIR
         DO J= 1,NVIR
            WAB(J,I)= WAB(J,I)-PAB(J,I)*(ENG(J+NOCC)+ENG(I+NOCC))
         ENDDO
      ENDDO
C
      CALL DGEMM('N','N',NBF,NOCC,NOCC,ONE,VECO2,NBF,
     *           PIJ,NOCC,ZER,PMX,NBF)
      CALL DGEMM('N','T',NBF,NBF,NOCC,ONE,PMX,NBF,
     *           VECO2,NBF,ZER,PMN,NBF)
      CALL DGEMM('N','N',NBF,NVIR,NVIR,-ONE,VECV,NBF,
     *           PAB,NVIR,ZER,PMX,NBF)
      CALL DGEMM('N','T',NBF,NBF,NVIR,ONE,PMX,NBF,
     *           VECV,NBF,ONE,PMN,NBF)
      CALL DAXPY(NOCC*NVAC,-HALF,WAI,1,XLAI,1)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT4A
      SUBROUTINE MPGT4A(NFT9,NFT15,NFT16,VECO,VECV,XAIBJC,XAIBJA,XKIBJC,
     *                  XKIBJA,DDIC,DDIA,PTR3C,PTR3A,LENDDI,LRECV,
     *                  NOCC,NVIR,NBF,NACORE,NOAC,MAXP,IPASS,MPROC,NUMI,
     *                  NUMV,NCYCV,MODV,NUMO,NCYCO,MODO,IRECO,IRECV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+0, ZER=0.0D+0)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION VECO(NBF,*),VECV(NBF,*),
     *          XAIBJC(NACORE,NOCC,NVIR),XAIBJA(NOAC,NOCC,NVIR),
     *          XKIBJC(NACORE,NOCC,NOCC),XKIBJA(NOAC,NOCC,NOCC),
     *          DDIC(NACORE*NOCC,MAXP,0:NPROC-1),
     *          DDIA(NOAC*NOCC,MAXP,0:NPROC-1),
     *          PTR3C(NACORE*NOCC,NBF),PTR3A(NOAC*NOCC,NBF),
     *          LENDDI(0:NPROC-1,3),LRECV(0:NPROC-1,2)
C
C     4TH QUARTER TRANSFORMATION, (AI|BJ) AND (KI|BJ)
C
      NOCAC=NOCC*NOAC
      NOCCO=NOCC*NACORE
C
      IF(GOPARR) THEN
      LDDIC=LENDDI(ME,1)
      LDDIA=LENDDI(ME,2)
      DO I= 1,NUMI
         IREC= (I-1)*NVIR+(IPASS-1)*NPROC+1
         DO IPROC = 0,MPROC-1
            READ(UNIT=NFT15,REC=IREC+IPROC)
     *          (DDIA(J,I,IPROC),J=1,NOCAC),(DDIC(J,I,IPROC),J=1,NOCCO)
         ENDDO
      ENDDO
C
      JPROC= ME
      DO IPROC= ME+1,NPROC-1
         JPROC= JPROC -1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDIC(1,1,IPROC),LDDIC*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3C(1,LRECV(JPROC,1)),LENDDI(JPROC,1)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
         CALL DDI_ISEND(DDIA(1,1,IPROC),LDDIA*8,IPROC,IREQ3)
         CALL DDI_IRECV(PTR3A(1,LRECV(JPROC,1)),LENDDI(JPROC,2)*8,JPROC,
     *                  IREQ4)
         CALL DDI_WAIT(IREQ3)
         CALL DDI_WAIT(IREQ4)
      ENDDO
      DO IPROC= 0,ME-1
         JPROC= JPROC-1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDIC(1,1,IPROC),LDDIC*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3C(1,LRECV(JPROC,1)),LENDDI(JPROC,1)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
         CALL DDI_ISEND(DDIA(1,1,IPROC),LDDIA*8,IPROC,IREQ3)
         CALL DDI_IRECV(PTR3A(1,LRECV(JPROC,1)),LENDDI(JPROC,2)*8,JPROC,
     *                  IREQ4)
         CALL DDI_WAIT(IREQ3)
         CALL DDI_WAIT(IREQ4)
      ENDDO
      CALL DCOPY(LDDIC,DDIC(1,1,ME),1,PTR3C(1,LRECV(ME,1)),1)
      CALL DCOPY(LDDIA,DDIA(1,1,ME),1,PTR3A(1,LRECV(ME,1)),1)
C
      IF(ME.GT.MPROC-1) RETURN
C
      NDIM= NUMO
      DO ICYC= 1,NCYCO
         IF(ICYC.EQ.NCYCO) NDIM= MODO
         IDIM= (ICYC-1)*NUMO+1
         IF(NACORE.NE.0)
     *   CALL DGEMM('N','N',NOCCO,NDIM,NBF,ONE,PTR3C,NOCCO,
     *              VECO(1,IDIM),NBF,ZER,XKIBJC(1,1,IDIM),NOCCO)
         CALL DGEMM('N','N',NOCAC,NDIM,NBF,ONE,PTR3A,NOCAC,
     *              VECO(1,IDIM),NBF,ZER,XKIBJA(1,1,IDIM),NOCAC)
         DO MOK= IDIM,IDIM+NDIM-1
            IRECO= IRECO+1
            WRITE(UNIT=NFT16,REC=IRECO)
     *           ((XKIBJC(J,I,MOK),J=1,NACORE),
     *            (XKIBJA(K,I,MOK),K=1,NOAC),I=NOAC+1,NOCC),
     *           ((XKIBJC(J,I,MOK),J=1,NACORE),
     *            (XKIBJA(K,I,MOK),K=1,NOAC),I=1,NOAC)
         ENDDO
      ENDDO
C
      NDIM= NUMV
      DO ICYC= 1,NCYCV
         IF(ICYC.EQ.NCYCV) NDIM= MODV
         IDIM= (ICYC-1)*NUMV+1
         IF(NACORE.NE.0)
     *   CALL DGEMM('N','N',NOCCO,NDIM,NBF,ONE,PTR3C,NOCCO,
     *              VECV(1,IDIM),NBF,ZER,XAIBJC(1,1,IDIM),NOCCO)
         CALL DGEMM('N','N',NOCAC,NDIM,NBF,ONE,PTR3A,NOCAC,
     *              VECV(1,IDIM),NBF,ZER,XAIBJA(1,1,IDIM),NOCAC)
         DO MOA= IDIM,IDIM+NDIM-1
            IRECV= IRECV+1
            WRITE(UNIT=NFT9, REC=IRECV)
     *           ((XAIBJC(J,I,MOA),J=1,NACORE),
     *            (XAIBJA(K,I,MOA),K=1,NOAC),I=NOAC+1,NOCC),
     *           ((XAIBJC(J,I,MOA),J=1,NACORE),
     *            (XAIBJA(K,I,MOA),K=1,NOAC),I=1,NOAC)
         ENDDO
      ENDDO
C
      ELSE
      DO I= 1,NBF
         IREC= (I-1)*NVIR+IPASS
         READ(UNIT=NFT15,REC=IREC)
     *       (PTR3A(J,I),J=1,NOCAC),(PTR3C(J,I),J=1,NOCCO)
      ENDDO
C
      NDIM= NUMO
      DO ICYC= 1,NCYCO
         IF(ICYC.EQ.NCYCO) NDIM= MODO
         IDIM= (ICYC-1)*NUMO+1
         IF(NACORE.NE.0)
     *   CALL DGEMM('N','N',NOCCO,NDIM,NBF,ONE,PTR3C,NOCCO,
     *              VECO(1,IDIM),NBF,ZER,XKIBJC(1,1,IDIM),NOCCO)
         CALL DGEMM('N','N',NOCAC,NDIM,NBF,ONE,PTR3A,NOCAC,
     *              VECO(1,IDIM),NBF,ZER,XKIBJA(1,1,IDIM),NOCAC)
         DO MOK= IDIM,IDIM+NDIM-1
            IRECO= IRECO+1
            WRITE(UNIT=NFT16,REC=IRECO)
     *           ((XKIBJC(J,I,MOK),J=1,NACORE),
     *            (XKIBJA(K,I,MOK),K=1,NOAC),I=NOAC+1,NOCC),
     *           ((XKIBJC(J,I,MOK),J=1,NACORE),
     *            (XKIBJA(K,I,MOK),K=1,NOAC),I=1,NOAC)
         ENDDO
      ENDDO
C
      NDIM= NUMV
      DO ICYC= 1,NCYCV
         IF(ICYC.EQ.NCYCV) NDIM= MODV
         IDIM= (ICYC-1)*NUMV+1
         IF(NACORE.NE.0)
     *   CALL DGEMM('N','N',NOCCO,NDIM,NBF,ONE,PTR3C,NOCCO,
     *              VECV(1,IDIM),NBF,ZER,XAIBJC(1,1,IDIM),NOCCO)
         CALL DGEMM('N','N',NOCAC,NDIM,NBF,ONE,PTR3A,NOCAC,
     *              VECV(1,IDIM),NBF,ZER,XAIBJA(1,1,IDIM),NOCAC)
         DO MOA= IDIM,IDIM+NDIM-1
            IRECV= IRECV+1
            WRITE(UNIT=NFT9, REC=IRECV)
     *           ((XAIBJC(J,I,MOA),J=1,NACORE),
     *            (XAIBJA(K,I,MOA),K=1,NOAC),I=NOAC+1,NOCC),
     *           ((XAIBJC(J,I,MOA),J=1,NACORE),
     *            (XAIBJA(K,I,MOA),K=1,NOAC),I=1,NOAC)
         ENDDO
      ENDDO
C
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGGEN
      SUBROUTINE MPGGEN(ENG,XAIBJA,XKIBJA,PAB,WAB,WAI,PIJ1,WIJ1,
     *                  TAB,DAIBJ,
     *                  NOCC,NVIR,NACORE,NOAC,NVAC,IPASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
      PARAMETER (MXRT=100)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
C
      DIMENSION ENG(*),XAIBJA(NOAC,NOCC,NVIR),XKIBJA(NOAC,NOCC,NOCC),
     *          PAB(*),WAB(*),WAI(*),PIJ1(*),WIJ1(*),
     *          TAB(NOAC,NOAC,NVAC),DAIBJ(NOAC,NOAC,NVAC)
C
      NOCAC=NOCC*NOAC
      NOAC2=NOAC*NOAC
C
      IF(GOPARR) THEN
         MOB= (IPASS-1)*NPROC+ME+1
         IF(MOB.GT.NVIR) RETURN
      ELSE
         MOB= IPASS
      ENDIF
      IF(MOB.GT.NVAC) GO TO 100
C
C     ----- CALCULATE MP2 ENERGY AND AMPLITUDE -----
C
      EB =-ENG(NOCC+MOB)
      DO MOA= 1,NVAC
         EAB =-ENG(NOCC+MOA)+EB
         DO MOI= 1,NOAC
            EABI = ENG(NACORE+MOI)+EAB
            DO MOJ= 1,NOAC
               EABIJ = ENG(NACORE+MOJ)+EABI
               ONEE= ONE/EABIJ
               TAB(MOJ,MOI,MOA)=(TWO*XAIBJA(MOJ,MOI,MOA)
     *                              -XAIBJA(MOI,MOJ,MOA))*ONEE
               DAIBJ(MOJ,MOI,MOA)= XAIBJA(MOJ,MOI,MOA)*ONEE
               E2 = E2+XAIBJA(MOJ,MOI,MOA)*TAB(MOJ,MOI,MOA)
            ENDDO
         ENDDO
      ENDDO
C
C     ----- CALCULATE PIJ AND WIJ[I] -----
C
      DO MOA= 1,NVAC
         CALL DGEMM('T','N',NOAC,NOAC,NOAC,ONE,TAB(1,1,MOA),NOAC,
     *              DAIBJ(1,1,MOA),NOAC,ONE,PIJ1,NOAC)
         CALL DGEMM('T','N',NOAC,NOCC,NOAC,-TWO,TAB(1,1,MOA),NOAC,
     *              XAIBJA(1,1,MOA),NOAC,ONE,WIJ1,NOAC)
      ENDDO
C
C     ----- CALCULATE WAI[I] -----
C
      CALL DGEMM('T','N',NOCC,NVAC,NOAC2,-FOUR,XKIBJA,NOCAC,
     *           TAB,NOAC2,ONE,WAI,NOCC)
C
C     ----- CALCULATE PAB AND WAB[I] -----
C
  100 CALL DGEMM('T','N',NVAC,NVAC,NOAC2,ONE,TAB,NOAC2,
     *           DAIBJ,NOAC2,ONE,PAB,NVIR)
      CALL DGEMM('T','N',NVAC,NVIR,NOAC2,-TWO,TAB,NOAC2,
     *           XAIBJA,NOCAC,ONE,WAB,NVAC)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGBT1
      SUBROUTINE MPGBT1(NFT20,VECV,TAB,TLB,TRECV,LENDDI,LRECV,
     *                  NVIR,NBF,NOAC,NVAC,MAXQ,IPASS,NUMS,MPROC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+0, ZER=0.0D+0)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION VECV(*),TAB(*),TLB(NOAC*NOAC,NBF),
     *          TRECV(NOAC*NOAC,MAXQ,0:NPROC-1),
     *          LENDDI(0:NPROC-1,3),LRECV(0:NPROC-1,2)
C
C     FIRST QUARTER BACK-TRANSFORMATION
C
      NOAC2=NOAC*NOAC
C
      IF(GOPARR) THEN
C
      MOB=(IPASS-1)*NPROC+ME+1
      IF(MOB.LE.NVAC) CALL DGEMM('N','N',NOAC2,NBF,NVAC,ONE,TAB,NOAC2,
     *                           VECV,NVIR,ZER,TLB,NOAC2)
      JPROC= ME
      DO IPROC= ME+1,NPROC-1
         JPROC= JPROC-1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(TLB(1,LRECV(IPROC,2)),LENDDI(IPROC,3)*8,IPROC,
     *                  IREQ1)
         CALL DDI_IRECV(TRECV(1,1,JPROC),LENDDI(ME,3)*8,JPROC,IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      DO IPROC= 0,ME-1
         JPROC= JPROC-1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(TLB(1,LRECV(IPROC,2)),LENDDI(IPROC,3)*8,IPROC,
     *                  IREQ1)
         CALL DDI_IRECV(TRECV(1,1,JPROC),LENDDI(ME,3)*8,JPROC,IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      CALL DCOPY(LENDDI(ME,3),TLB(1,LRECV(ME,2)),1,TRECV(1,1,ME),1)
C
      DO I= 1,NUMS
         IREC=(I-1)*NBF+(IPASS-1)*NPROC
         DO J= 1,MPROC
            WRITE(UNIT=NFT20,REC=IREC+J) (TRECV(K,I,J-1),K=1,NOAC2)
         ENDDO
      ENDDO
C
      ELSE
C
      IF(IPASS.GT.NVAC) RETURN
      CALL DGEMM('N','N',NOAC2,NBF,NVAC,ONE,TAB,NOAC2,
     *           VECV,NVIR,ZER,TLB,NOAC2)
C
      DO I= 1,NBF
        IREC=(I-1)*NBF+IPASS
        WRITE(UNIT=NFT20,REC=IREC) (TLB(J,I),J=1,NOAC2)
      ENDDO
C
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGQTC
      SUBROUTINE MPGQTC(NFT15,NFT19,NFT20,VECO,VECO2,VECO3,VECV,VECV2,
     *                  XLAI,XLNI1,XLNI2,GHONDO,XINTS,AOINTS,
     *                  PTR1,PTR2,PTR3,TLB,TLS,PMN,
     *                  NOCC,NVIR,NBF,NOAC,NVAC,NACORE,NANGM,CUTOFF,
     *                  NLOOP,ILOOP,NA3,NS1,NBT2,DIVANG,MXANGM,
     *                  LBLINT,NUMINT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL DIVANG,SCHSKP
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (ONE=1.0D+00)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C     COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
C
      DIMENSION VECO(*),VECO2(*),VECO3(*),VECV(*),VECV2(*),XLAI(NOCC,*),
     *          XLNI1(NOCC,*),XLNI2(NOCC,*),GHONDO(*),XINTS(*),
     *          AOINTS(*),PTR1(*),PTR2(*),PTR3(*),TLB(*),TLS(*),PMN(*)
      DIMENSION NLOOP(0:NPROC-1,4),ILOOP(NSHELL,0:NPROC-1,2)
      DIMENSION LBLINT(NBF*(NANGM**3)),NUMINT(NANGM**3)
C
C     1ST, 2ND, 3RD QUARTER TRANSFORMATIONS
C
      ISTART= 1
      IEND  = NSHELL
      ICOUNT= 0
      IRECO = 0
      IRECV = 0
      NOCC3 = (NOCC*NOCC+NOCC)/2
      IF(GOPARR) IEND= NLOOP(ME,3)
C
      CALL VCLR(XLNI1,1,NOCC*NBF)
      CALL VCLR(XLNI2,1,NOCC*NBF)
C
C     ----- START I SHELL -----
C
      DO 500 II= ISTART,IEND
C
C     ----- GO PARALLEL! -----
C
         ISH = II
         IF(GOPARR) ISH = ILOOP(II,ME,2)
C
         MINI= KMIN(ISH)
         MAXI= KMAX(ISH)
         LOCI= KLOC(ISH)-MINI
         LOCIE=LOCI+MAXI
         MINI2= MINI
         MAXI2= MAXI
         ITMP= MAXI2-MINI2+1
         IDIV= 0
  110    IF(DIVANG.AND.ITMP.GT.MXANGM) THEN
            MINI2= MXANGM*IDIV+KMIN(ISH)
            MAXI2= MXANGM+MINI2-1
            IF(MAXI2.GT.KMAX(ISH)) MAXI2=KMAX(ISH)
            IDIV= IDIV+1
         ENDIF
         NUMI= MAXI2-MINI2+1
         CALL VCLR(PTR2,1,NOCC3*NBF*NUMI)
         JSTART= 0
         JCOUNT= 0
         CALL MPGBT2(NFT20,VECV2,TLB,TLS,NVIR,NBF,NOAC,NVAC,NUMI,NBT2,
     *               ICOUNT,LOCIE)
C
C     ----- START J SHELL -----
C
      DO 400 JSH= 1,NSHELL
         IJIJ= IA(ISH)+JSH
         IF(JSH.GT.ISH) IJIJ= IA(JSH)+ISH
         MINJ= KMIN(JSH)
         MAXJ= KMAX(JSH)
         NUMJ= MAXJ-MINJ+1
         IF(JCOUNT.EQ.0) CALL VCLR(PTR1,1,NOCC*NS1*NBF*NUMI)
C
C     ----- START K SHELL -----
C
      DO 300 KSH= 1,NSHELL
         MINK= KMIN(KSH)
         MAXK= KMAX(KSH)
         NUMK= MAXK-MINK+1
         LOCK= KLOC(KSH)-MINK
         DO ITMP= 1,NUMI*NUMJ*NUMK
            NUMINT(ITMP)=0
         ENDDO
C
C     ----- START L SHELL -----
C
      DO 200 LSH= 1,KSH
         MINL= KMIN(LSH)
         MAXL= KMAX(LSH)
         LOCL= KLOC(LSH)-MINL
         KLKL= IA(KSH)+LSH
         TEST= XINTS(IJIJ)*XINTS(KLKL)
         SCHSKP= TEST.LT.CUTOFF
C
         IF(.NOT.SCHSKP) THEN
            CALL SHELLQUARTMP2(GHONDO)
C
C     ----- SAVE AO INTEGRALS -----
C
            CALL SAVEAO2(AOINTS,GHONDO,CUTOFF,
     *                   LBLINT,NUMINT,KSH,LSH,LOCL,MINI,NBF,
     *                   MINI2,MAXI2,MINJ,MAXJ,MINK,MAXK,MINL,MAXL)
         ENDIF
 200  CONTINUE
C
C     ----- END L SHELL -----
C
      CALL MPGT1B(VECO,AOINTS,PTR1,LBLINT,NUMINT,NOCC,NBF,NS1,
     *            NUMI,NUMJ,LOCK,MINK,MAXK,JCOUNT)
 300  CONTINUE
C
C     ----- END K SHELL -----
C
      JCOUNT= JCOUNT+NUMJ
      IF(JSH.NE.NSHELL) JNEXT= JCOUNT+KMAX(JSH+1)-KMIN(JSH+1)+1
      IF(JSH.NE.NSHELL.AND.JNEXT.LE.NS1) GO TO 400
C
      CALL MPGT2B(VECO2,PTR1,PTR2,NOCC,NBF,NOAC,NACORE,NOCC3,
     *            NS1,NUMI,JSTART,JCOUNT,CUTOFF)
      CALL MPGLNI(PTR1,TLS,PMN,XLNI1,NOCC,NBF,NOAC,NS1,
     *            NUMI,LOCI,MINI2,JSTART,JCOUNT)
C
      JSTART= JSTART+JCOUNT
      JCOUNT= 0
 400  CONTINUE
C
C     ----- END J SHELL -----
C
      CALL MPGT3B(NFT15,NFT19,VECO3,VECV,PTR2,PTR3,NOCC,NVIR,NBF,NOCC3,
     *            NUMI,NA3,IRECO,IRECV)
C
      ICOUNT= ICOUNT+NUMI
      IF(DIVANG.AND.MAXI2.NE.KMAX(ISH)) GO TO 110
C
 500  CONTINUE
C
C     ----- END I SHELL -----
C
      DO LN= 1,NBF
         DO LI= 1,NACORE
            XLNI2(LI,LN)= XLNI1(LI+NOAC,LN)
         ENDDO
         DO LI= 1,NOAC
            XLNI2(LI+NACORE,LN)= XLNI1(LI,LN)
         ENDDO
      ENDDO
      IF(GOPARR) CALL DDI_GSUMF(1001,XLNI2,NOCC*NBF)
C
      CALL DGEMM('N','N',NOCC,NVIR,NBF,ONE,XLNI2,NOCC,
     *           VECV,NBF,ONE,XLAI,NOCC)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT1B
      SUBROUTINE MPGT1B(VECO,AOINTS,PTR1,LBLINT,NUMINT,NOCC,NBF,NS1,
     *                  NUMI,NUMJ,LOCK,MINK,MAXK,JCOUNT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VECO(NOCC,*),AOINTS(NBF,*),PTR1(NOCC,NS1,NBF,*),
     *          LBLINT(NBF,*),NUMINT(*)
C
C     1ST QUARTER TRANSFORMATION
C
      IJK= 0
      DO I= 1,NUMI
         DO J= 1,NUMJ
            JJ= JCOUNT+J
            DO K= MINK,MAXK
               KK= LOCK+K
               IJK= IJK+1
               DO L=1,NUMINT(IJK)
                  LL = LBLINT(L,IJK)
                  DUM= AOINTS(L,IJK)
                  DO MOI= 1,NOCC
                     PTR1(MOI,JJ,KK,I)=PTR1(MOI,JJ,KK,I)
     *                                +DUM*VECO(MOI,LL)
                  ENDDO
                  DO MOI= 1,NOCC
                     PTR1(MOI,JJ,LL,I)=PTR1(MOI,JJ,LL,I)
     *                                +DUM*VECO(MOI,KK)
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT2B
      SUBROUTINE MPGT2B(VECO,PTR1,PTR2,NOCC,NBF,NOAC,NACORE,NOCC3,
     *                  NS1,NUMI,JSTART,JCOUNT,CUTOFF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VECO(NOCC,*),PTR1(NOCC,NS1,NBF,*),PTR2(NOCC3,NBF,*)
C
C     2ND QUARTER TRANSFORMATION
C
      DO 190 I= 1,NUMI
         DO 180 J= 1,JCOUNT
            JJ= JSTART+J
            DO 170 K= 1,NBF
               MOIJ= 0
               DO 130 MOI= 1,NACORE
                  DUM= PTR1(MOI+NOAC,J,K,I)
                  IF(ABS(DUM).LT.CUTOFF) GO TO 120
                  DO 110 MOJ= 1,MOI
                     PTR2(MOJ+MOIJ,JJ,I)= PTR2(MOJ+MOIJ,JJ,I)
     *                                   +DUM*VECO(MOJ,K)
  110             CONTINUE
  120             MOIJ= MOIJ+MOI
  130          CONTINUE
               DO 160 MOI= 1,NOAC
                  DUM= PTR1(MOI,J,K,I)
                  IF(ABS(DUM).LT.CUTOFF) GO TO 150
                  DO 140 MOJ= 1,MOI+NACORE
                     PTR2(MOJ+MOIJ,JJ,I)= PTR2(MOJ+MOIJ,JJ,I)
     *                                   +DUM*VECO(MOJ,K)
  140             CONTINUE
  150             MOIJ= MOIJ+MOI+NACORE
  160          CONTINUE
  170       CONTINUE
  180    CONTINUE
  190 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT3B
      SUBROUTINE MPGT3B(NFT15,NFT19,VECO,VECV,PTR2,PTR3,NOCC,NVIR,NBF,
     *                  NOCC3,NUMI,NA3,IRECO,IRECV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
C
      DIMENSION VECO(NBF,*),VECV(NBF,*),PTR2(NOCC3*NBF,*),PTR3(NOCC3,*)
C
C     3RD QUARTER TRANSFORMATION
C
      NPASSO= NOCC/NA3
      NPASSV= NVIR/NA3
      MODO  = MOD(NOCC,NA3)
      MODV  = MOD(NVIR,NA3)
      IF(MODO.NE.0) THEN
         NPASSO= NPASSO+1
      ELSE
         MODO= NA3
      ENDIF
      IF(MODV.NE.0) THEN
         NPASSV= NPASSV+1
      ELSE
         MODV= NA3
      ENDIF
C
      DO I= 1,NUMI
         MK3= NA3
         DO IPASS= 1,NPASSO
            IF(IPASS.EQ.NPASSO) MK3= MODO
            MOK= (IPASS-1)*NA3+1
            CALL DGEMM('N','N',NOCC3,MK3,NBF,ONE,PTR2(1,I),NOCC3,
     *                 VECO(1,MOK),NBF,ZER,PTR3,NOCC3)
            DO IWRITE= 1,MK3
               IRECO= IRECO+1
               WRITE(UNIT=NFT15,REC=IRECO)
     *              (PTR3(MOIJ,IWRITE),MOIJ=1,NOCC3)
            ENDDO
         ENDDO
      ENDDO
C
      DO I= 1,NUMI
         MA3= NA3
         DO IPASS= 1,NPASSV
            IF(IPASS.EQ.NPASSV) MA3= MODV
            MOA= (IPASS-1)*NA3+1
            CALL DGEMM('N','N',NOCC3,MA3,NBF,ONE,PTR2(1,I),NOCC3,
     *                 VECV(1,MOA),NBF,ZER,PTR3,NOCC3)
            DO IWRITE= 1,MA3
               IRECV= IRECV+1
               WRITE(UNIT=NFT19,REC=IRECV)
     *              (PTR3(MOIJ,IWRITE),MOIJ=1,NOCC3)
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK SAVEAO2
      SUBROUTINE SAVEAO2(AOINTS,GHONDO,CUTOFF,
     *                   LBLINT,NUMINT,KSH,LSH,LOCL,MINI,NBF,
     *                   MINI2,MAXI2,MINJ,MAXJ,MINK,MAXK,MINL,MAXL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL KANDL
C
      PARAMETER (HALF=0.5D+00)
C
      COMMON /ERIOUT/ INU,JNU,KNU,LNU,LSTRI,LSTRJ,LSTRK,LSTRL
C
      DIMENSION AOINTS(NBF,*),GHONDO(*)
      DIMENSION LBLINT(NBF,*),NUMINT(*)
C
C     SAVE AO INTEGRALS
C
      KANDL= KSH.EQ.LSH
      IJK  = 0
      LMAX = MAXL
C
      DO 240 I= MINI2,MAXI2
         I_INDEX = (I-MINI)*LSTRI + 1
         DO 230 J= MINJ,MAXJ
            IJ_INDEX = (J-MINJ)*LSTRJ + I_INDEX
            DO 220 K= MINK,MAXK
               IJK_INDEX = (K-MINK)*LSTRK + IJ_INDEX
               IJK= IJK+1
               IF(KANDL) LMAX= K
               DO 210 L= MINL,LMAX
                  LL= LOCL+L
                  IJKL_INDEX = (L-MINL)*LSTRL + IJK_INDEX
C
                  VAL = GHONDO(IJKL_INDEX)
                  IF(ABS(VAL).LT.CUTOFF) GO TO 210
                  IF(KANDL.AND.K.EQ.L) VAL= VAL*HALF
                  NUMINT(IJK)= NUMINT(IJK)+1
                  AOINTS(NUMINT(IJK),IJK)= VAL
                  LBLINT(NUMINT(IJK),IJK)= LL
  210          CONTINUE
  220       CONTINUE
  230    CONTINUE
  240 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGBT2
      SUBROUTINE MPGBT2(NFT20,VECV,TLB,TLS,NVIR,NBF,NOAC,NVAC,NUMI,NBT2,
     *                  ICOUNT,LOCIE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ONE=1.0D+0, ZER=0.0D+0)
C
      DIMENSION VECV(NVIR,*),TLB(NOAC*NOAC,NVAC),TLS(NOAC*NOAC,NBF,*)
C
C     2ND BACK-QUARTER TRANSFORMATION
C
      NOAC2= NOAC*NOAC
      NCYCLE= NBF/NBT2
      MODBT2= MOD(NBF,NBT2)
      IF(MODBT2.NE.0) THEN
         NCYCLE= NCYCLE+1
      ELSE
         MODBT2= NBT2
      ENDIF
C
      DO II= 1,NUMI
         IREC= (ICOUNT+II-1)*NBF
         DO MOB= 1,NVAC
            READ(UNIT=NFT20,REC=IREC+MOB) (TLB(MOIJ,MOB),MOIJ=1,NOAC2)
         ENDDO
         MBT2= NBT2
         DO ICYCLE= 1,NCYCLE
            MS= (ICYCLE-1)*NBT2+1
            IF(ICYCLE.EQ.NCYCLE) MBT2= MODBT2
            CALL DGEMM('N','N',NOAC2,MBT2,NVAC,ONE,TLB,NOAC2,
     *                 VECV(1,MS),NVIR,ZER,TLS(1,MS,II),NOAC2)
C           IF(MS.LT.LOCIE) THEN
               DO IS= MS,MS+MBT2-1
                  IF(IS.LE.LOCIE) THEN
                  IREC= IREC+1
                  WRITE(UNIT=NFT20,REC=IREC)
     *                 (TLS(MOIJ,IS,II),MOIJ=1,NOAC2)
                  ENDIF
               ENDDO
C           ENDIF
         ENDDO
      ENDDO
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGLNI
      SUBROUTINE MPGLNI(PTR1,TLS,PMN,XLNI1,NOCC,NBF,NOAC,NS1,
     *                  NUMI,LOCI,MINI2,JSTART,JCOUNT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
C
      DIMENSION PTR1(NOCC*NS1,NBF,*),TLS(NOAC*NOAC,NBF,*),PMN(NBF,*),
     *          XLNI1(NOCC,*)
C
C     CALCULATE MP2 LAGRANGIAN
C
      II= LOCI+MINI2
      CALL DGEMV('N',NOCC*JCOUNT,NBF*NUMI,-TWO,PTR1,NOCC*NS1,
     *           PMN(1,II),1,ONE,XLNI1(1,JSTART+1),1)
C
      DO I= 1,NUMI
         DO K= 1,NBF
            CALL DGEMV('N',NOCC,JCOUNT,FOUR,PTR1(1,K,I),NOCC,
     *                 PMN(JSTART+1,II),1,ONE,XLNI1(1,K),1)
            CALL DGEMM('T','N',NOAC,JCOUNT,NOAC,-TWO,TLS(1,K,I),NOAC,
     *                 PTR1(1,K,I),NOCC,ONE,XLNI1(1,JSTART+1),NOCC)
C    *                 PTR1(1,K,I),NOCC,ONE,XLNI2(1,JSTART+1),NOCC)
         ENDDO
         II= II+1
      ENDDO
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGQTD
      SUBROUTINE MPGQTD(NFT9,NFT15,NFT19,VECO,VECV,ENG,PIJ,PAB,PTR3,
     *                  DDI,PTRO4,PTRV4,AIBJ,WIJ,WIJ1,
     *                  WIJ2,WORK,ENGAI,NOCC,NVIR,NBF,NOCC3,NLOOP,
     *                  MAXQ,NWORK,LENDDI,LRECV,NPROCS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00, MXAO=8192)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IJPAIR/ IA(MXAO)
C
      DIMENSION VECO(*),VECV(*),ENG(*),PIJ(NOCC,*),PAB(NVIR,*),PTR3(*),
     *          DDI(*),PTRO4(NOCC3,*),PTRV4(NOCC3,*),AIBJ(NOCC,NOCC,*),
     *          WIJ(*),WIJ1(*),WIJ2(NOCC,*),WORK(NOCC,NOCC,*),
     *          ENGAI(NOCC,*),NLOOP(0:NPROC-1,4)
      DIMENSION LENDDI(0:NPROCS-1),LRECV(0:NPROCS-1)
C
      NOCC2= NOCC*NOCC
      CALL VCLR(WIJ1,1,NOCC3)
      CALL VCLR(WIJ2,1,NOCC2)
C
C     ----- CALCULATE SIZE AND POINTER FOR DDI_ISEND AND DDI_IRECV -----
C
      IF(GOPARR) THEN
         NUMI= NLOOP(ME,4)
         NSUM= 1
         DO IPROC= 0,NPROC-1
            LENDDI(IPROC)= NLOOP(IPROC,4)*NOCC3
            LRECV(IPROC) = NSUM
            NSUM= NSUM+NLOOP(IPROC,4)
         ENDDO
      ENDIF
      DO MOA= 1,NVIR
         ENGA= ENG(MOA+NOCC)
         DO MOI= 1,NOCC
            ENGAI(MOI,MOA)= ONE/(ENGA-ENG(MOI))
         ENDDO
      ENDDO
C
C     ----- 4TH QUARTER TRANSFORMATION OF (IJ|KL) -----
C
      NPASS= NOCC/NPROC
      IF(MOD(NOCC,NPROC).NE.0) NPASS= NPASS+1
C
      DO 110 IPASS= 1,NPASS
         MPROC= NPROC
         IF(IPASS.EQ.NPASS) THEN
            MPROC= MOD(NOCC,NPROC)
            IF(MPROC.EQ.0) MPROC= NPROC
         ENDIF
C
         CALL MPGT4B(NFT15,VECO,PTR3,DDI,PTRO4,LENDDI,LRECV,NOCC,NOCC3,
     *               NBF,IPASS,NUMI,MPROC,MAXQ)
C
C     ----- CALCULATE A PART OF WIJ[III] -----
C
         MOK= IPASS
         IF(GOPARR) MOK=(IPASS-1)*NPROC+ME+1
         IF(MOK.LE.NOCC) THEN
            CALL DGEMV('N',NOCC3,NOCC,FOUR,PTRO4,NOCC3,
     *                 PIJ(1,MOK),1,ONE,WIJ1,1)
            DO MOL= 1,NOCC
               CALL DSPMV('U',NOCC,-TWO,PTRO4(1,MOL),PIJ(1,MOK),1,ONE,
     *                    WIJ2(1,MOL),1)
            ENDDO
         ENDIF
  110 CONTINUE
C
C     ----- 4TH QUARTER TRANSFORMATION OF (IJ|AB) -----
C
      NPASS= NVIR/NPROC
      IF(MOD(NVIR,NPROC).NE.0) NPASS= NPASS+1
      DO 120 IPASS= 1,NPASS
         MPROC= NPROC
         IF(IPASS.EQ.NPASS) THEN
            MPROC= MOD(NVIR,NPROC)
            IF(MPROC.EQ.0) MPROC= NPROC
         ENDIF
C
         CALL MPGT4C(NFT19,VECV,PTR3,DDI,PTRV4,LENDDI,LRECV,NVIR,NOCC3,
     *               NBF,IPASS,NUMI,MPROC,MAXQ)
C
         MOA= IPASS
         IF(GOPARR) MOA=(IPASS-1)*NPROC+ME+1
         IF(MOA.LE.NVIR) THEN
C
C     ----- READ (AI|BJ) -----
C
            IREC=(IPASS-1)*NVIR
            DO I= 1,NVIR
               READ(UNIT=NFT9,REC=IREC+I)
     *             ((AIBJ(K,J,I),K=1,NOCC),J=1,NOCC)
            ENDDO
C
C     ----- CALCULATE AND SAVE 4(AI|BJ)-(AB|IJ)-(AJ|BI) -----
C
            NCYCLE= NVIR/NWORK
            MODW= MOD(NVIR,NWORK)
            IF(MODW.NE.0) THEN
               NCYCLE= NCYCLE+1
            ELSE
               MODW= NWORK
            ENDIF
            MWORK= NWORK
            DO ICYCLE= 1,NCYCLE
               IF(ICYCLE.EQ.NCYCLE) MWORK= MODW
               NOB= (ICYCLE-1)*NWORK
               DO IB= 1,MWORK
                  MOB= NOB+IB
                  DO I= 1,NOCC
                     DO J= 1,NOCC
                        IJ=IA(I)+J
                        IF(J.GT.I) IJ= IA(J)+I
                        WORK(J,I,IB)=(FOUR*AIBJ(J,I,MOB)-AIBJ(I,J,MOB)
     *                               -PTRV4(IJ,MOB))*ENGAI(J,MOA)
                     ENDDO
                  ENDDO
                  IF(MOB.EQ.MOA) THEN
                     DO I=1,NOCC
                        WORK(I,I,IB)= WORK(I,I,IB)+ONE
                     ENDDO
                  ENDIF
               ENDDO
               DO IB= 1,MWORK
                  WRITE(UNIT=NFT9,REC=IREC+IB)
     *                 ((WORK(J,I,IB),J=1,NOCC),I=1,NOCC)
               ENDDO
C
               IREC= IREC+MWORK
            ENDDO
C
C     ----- CALCULATE A PART OF WIJ[III] -----
C
            CALL DGEMV('N',NOCC2,NVIR,TWO,AIBJ,NOCC2,
     *                 PAB(1,MOA),1,ONE,WIJ2,1)
            CALL DGEMV('N',NOCC3,NVIR,-FOUR,PTRV4,NOCC3,
     *                 PAB(1,MOA),1,ONE,WIJ1,1)
         ENDIF
  120 CONTINUE
C
      II= 0
      DO I= 1,NOCC
         DO J= 1,I-1
            II= II+1
            WIJ2(I,J)= WIJ2(I,J)+WIJ1(II)
            WIJ2(J,I)= WIJ2(J,I)+WIJ1(II)
         ENDDO
         II= II+1
         WIJ2(I,I)= WIJ2(I,I)+WIJ1(II)
      ENDDO
      IF(GOPARR) CALL DDI_GSUMF(1002,WIJ2,NOCC2)
      CALL DAXPY(NOCC2,ONE,WIJ2,1,WIJ,1)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT4B
      SUBROUTINE MPGT4B(NFT15,VECO,PTR3,DDI,PTR4,LENDDI,LRECV,NOCC,
     *                  NOCC3,NBF,IPASS,NUMI,MPROC,MAXQ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, ZER=0.0D+00)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION VECO(*),PTR3(NOCC3,*),DDI(NOCC3,MAXQ,0:NPROC-1),PTR4(*),
     *          LENDDI(0:NPROC-1),LRECV(0:NPROC-1)
C
C     4TH QUARTER TRANSFORMATION OF (IJ|KL)
C
      IF(GOPARR) THEN
      LDDI= LENDDI(ME)
      DO I= 1,NUMI
         IREC= (I-1)*NOCC+(IPASS-1)*NPROC+1
         DO IPROC= 0,MPROC-1
            READ(UNIT=NFT15,REC=IREC+IPROC) (DDI(J,I,IPROC),J=1,NOCC3)
         ENDDO
      ENDDO
C
      JPROC= ME
      DO IPROC= ME+1,NPROC-1
         JPROC= JPROC -1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDI(1,1,IPROC),LDDI*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3(1,LRECV(JPROC)),LENDDI(JPROC)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      DO IPROC= 0, ME-1
         JPROC= JPROC -1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDI(1,1,IPROC),LDDI*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3(1,LRECV(JPROC)),LENDDI(JPROC)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      CALL DCOPY(LDDI,DDI(1,1,ME),1,PTR3(1,LRECV(ME)),1)
C
      IF(ME.GT.MPROC-1) RETURN
C
      CALL DGEMM('N','N',NOCC3,NOCC,NBF,ONE,PTR3,NOCC3,
     *           VECO,NBF,ZER,PTR4,NOCC3)
C
      ELSE
      DO I= 1,NBF
         IREC= (I-1)*NOCC+IPASS
         READ(UNIT=NFT15,REC=IREC) (PTR3(J,I),J=1,NOCC3)
      ENDDO
C
      CALL DGEMM('N','N',NOCC3,NOCC,NBF,ONE,PTR3,NOCC3,
     *           VECO,NBF,ZER,PTR4,NOCC3)
C
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGT4C
      SUBROUTINE MPGT4C(NFT19,VECV,PTR3,DDI,PTR4,LENDDI,LRECV,NVIR,
     *                  NOCC3,NBF,IPASS,NUMI,MPROC,MAXQ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, ZER=0.0D+00)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION VECV(*),PTR3(NOCC3,*),DDI(NOCC3,MAXQ,0:NPROC-1),PTR4(*),
     *          LENDDI(0:NPROC-1),LRECV(0:NPROC-1)
C
C     4TH QUARTER TRANSFORMATION OF (IJ|AB)
C
      IF(GOPARR) THEN
      LDDI= LENDDI(ME)
      DO I= 1,NUMI
         IREC= (I-1)*NVIR+(IPASS-1)*NPROC+1
         DO IPROC= 0,MPROC-1
            READ(UNIT=NFT19,REC=IREC+IPROC) (DDI(J,I,IPROC),J=1,NOCC3)
         ENDDO
      ENDDO
C
      JPROC= ME
      DO IPROC= ME+1,NPROC-1
         JPROC= JPROC -1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDI(1,1,IPROC),LDDI*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3(1,LRECV(JPROC)),LENDDI(JPROC)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      DO IPROC= 0, ME-1
         JPROC= JPROC -1
         IF(JPROC.EQ.-1) JPROC= NPROC-1
         CALL DDI_ISEND(DDI(1,1,IPROC),LDDI*8,IPROC,IREQ1)
         CALL DDI_IRECV(PTR3(1,LRECV(JPROC)),LENDDI(JPROC)*8,JPROC,
     *                  IREQ2)
         CALL DDI_WAIT(IREQ1)
         CALL DDI_WAIT(IREQ2)
      ENDDO
      CALL DCOPY(LDDI,DDI(1,1,ME),1,PTR3(1,LRECV(ME)),1)
C
      IF(ME.GT.MPROC-1) RETURN
C
      CALL DGEMM('N','N',NOCC3,NVIR,NBF,ONE,PTR3,NOCC3,
     *           VECV,NBF,ZER,PTR4,NOCC3)
C
      ELSE
      DO I= 1,NBF
         IREC= (I-1)*NVIR+IPASS
         READ(UNIT=NFT19,REC=IREC) (PTR3(J,I),J=1,NOCC3)
      ENDDO
C
      CALL DGEMM('N','N',NOCC3,NVIR,NBF,ONE,PTR3,NOCC3,
     *           VECV,NBF,ZER,PTR4,NOCC3)
C
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGCPHF
      SUBROUTINE MPGCPHF(NFT9,ENG,PAI,XLAI,AIBJ,GMRES1,
     *                   GMRES2,GMVEC,GMTRN,GMWRK1,GMWRK2,
     *                   GMFAC,NOCC,NVIR,ITERMX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, ZER=0.0D+00)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION ENG(*),PAI(NOCC,*),XLAI(NOCC,*),AIBJ(NOCC*NOCC,*),
     *          GMRES1(ITERMX+1,*),GMRES2(ITERMX+1,*),
     *          GMVEC(NOCC,NVIR,0:ITERMX+1),GMTRN(ITERMX+1,*),
     *          GMWRK1(ITERMX+1,*),GMWRK2(ITERMX+1,*),GMFAC(*)
C
C     SOLVE THE CPHF EQUATIONS OF MP2 USING GMRES METHOD
C
      THRE= 1.D-10
      CALL VCLR(GMRES1,1,ITERMX*(ITERMX+1))
      CALL VCLR(GMRES2,1,ITERMX*(ITERMX+1))
      CALL VCLR(GMVEC(1,1,0),1,NOCC*NVIR*(ITERMX+2))
      CALL VCLR(GMTRN,1,(ITERMX+1)*(ITERMX+1))
      CALL VCLR(GMWRK1,1,(ITERMX+1)*(ITERMX+1))
      DO II= 1,ITERMX+1
         GMTRN(II,II)= ONE
         GMWRK1(II,II)= ONE
      ENDDO
C
C     ----- CALCULATE INITIAL VECTOR -----
C
      DO MOA= 1,NVIR
         ENGA= ENG(MOA+NOCC)
         DO MOI= 1,NOCC
            XLAI(MOI,MOA)= XLAI(MOI,MOA)/(ENGA-ENG(MOI))
         ENDDO
      ENDDO
C
      NPASS= NVIR/NPROC
      IF(MOD(NVIR,NPROC).NE.0) NPASS= NPASS+1
      CALL DCOPY(NOCC*NVIR,XLAI,1,GMVEC(1,1,0),1)
C
      DO IPASS= 1,NPASS
         IF(GOPARR) THEN
            MOB=(IPASS-1)*NPROC+ME+1
            IF(MOB.GT.NVIR) GO TO 300
         ELSE
            MOB= IPASS
         ENDIF
         ICOUNT=(IPASS-1)*NVIR
         CALL DCOPY(NOCC,XLAI(1,MOB),1,GMVEC(1,MOB,1),1)
         DO MOA= 1,NVIR
           READ(UNIT=NFT9,REC=ICOUNT+MOA) (AIBJ(I,MOA),I=1,NOCC*NOCC)
         ENDDO
         CALL DGEMV('N',NOCC,NOCC*NVIR,-ONE,AIBJ,NOCC,
     *              GMVEC(1,1,0),1,ONE,GMVEC(1,MOB,1),1)
      ENDDO
  300 IF(GOPARR) CALL DDI_GSUMF(1001,GMVEC(1,1,1),NOCC*NVIR)
      GAMMA= DDOT(NOCC*NVIR,GMVEC(1,1,1),1,GMVEC(1,1,1),1)
      GAMMA= SQRT(GAMMA)
      CALL DSCAL(NOCC*NVIR,ONE/GAMMA,GMVEC(1,1,1),1)
C
C     ----- START ITERATIONS -----
C
      DO ICYCLE= 1,ITERMX
         DO IPASS= 1,NPASS
            IF(GOPARR) THEN
               MOB=(IPASS-1)*NPROC+ME+1
               IF(MOB.GT.NVIR) GO TO 200
            ELSE
               MOB= IPASS
            ENDIF
C
            ICOUNT=(IPASS-1)*NVIR
            DO MOA= 1,NVIR
              READ(UNIT=NFT9,REC=ICOUNT+MOA) (AIBJ(I,MOA),I=1,NOCC*NOCC)
            ENDDO
C
            CALL DGEMV('N',NOCC,NOCC*NVIR,ONE,AIBJ,NOCC,
     *                 GMVEC(1,1,ICYCLE),1,ZER,GMVEC(1,MOB,ICYCLE+1),1)
         ENDDO
  200    IF(GOPARR) CALL DDI_GSUMF(1002,GMVEC(1,1,ICYCLE+1),NOCC*NVIR)
C
         DO II= 1,ICYCLE
            GMRES1(II,ICYCLE)= DDOT(NOCC*NVIR,GMVEC(1,1,II),1,
     *                              GMVEC(1,1,ICYCLE+1),1)
         ENDDO
C
         DO II= 1,ICYCLE
            CALL DAXPY(NOCC*NVIR,-GMRES1(II,ICYCLE),GMVEC(1,1,II),1,
     *                 GMVEC(1,1,ICYCLE+1),1)
         ENDDO
C
         FAC=DDOT(NOCC*NVIR,GMVEC(1,1,ICYCLE+1),1,GMVEC(1,1,ICYCLE+1),1)
         FAC= SQRT(FAC)
         GMRES1(ICYCLE+1,ICYCLE)= FAC
         FAC= ONE/FAC
         CALL DSCAL(NOCC*NVIR,FAC,GMVEC(1,1,ICYCLE+1),1)
C
         CALL DGEMM('N','N',ICYCLE+1,ICYCLE,ICYCLE+1,ONE,GMTRN,
     *              ITERMX+1,GMRES1,ITERMX+1,ZER,GMRES2,ITERMX+1)
         FAC= ONE/(GMRES2(ICYCLE,ICYCLE)**2+GMRES2(ICYCLE+1,ICYCLE)**2)
         FAC= SQRT(FAC)
         WORK1= GMRES2(ICYCLE,ICYCLE)*FAC
         WORK2= GMRES2(ICYCLE+1,ICYCLE)*FAC
         GMWRK1(ICYCLE,ICYCLE)= WORK1
         GMWRK1(ICYCLE+1,ICYCLE)=-WORK2
         GMWRK1(ICYCLE,ICYCLE+1)= WORK2
         GMWRK1(ICYCLE+1,ICYCLE+1)= WORK1
         CALL DGEMM('N','N',ICYCLE+1,ICYCLE+1,ICYCLE+1,ONE,GMWRK1,
     *              ITERMX+1,GMTRN,ITERMX+1,ZER,GMWRK2,ITERMX+1)
         DO II= 1,ICYCLE+1
            DO JJ= 1,ICYCLE+1
               GMTRN(JJ,II) =GMWRK2(JJ,II)
            ENDDO
         ENDDO
         GMWRK1(ICYCLE,ICYCLE)= ONE
         GMWRK1(ICYCLE+1,ICYCLE)= ZER
         GMWRK1(ICYCLE,ICYCLE+1)= ZER
         GMWRK1(ICYCLE+1,ICYCLE+1)= ONE
C
         NCYCLE= ICYCLE
         IF(MASWRK) WRITE(IW,9100) ICYCLE,GMTRN(ICYCLE+1,1)*GAMMA
         IF(ABS(GMTRN(ICYCLE+1,1)*GAMMA).LT.THRE) GO TO 100
      ENDDO
C
  100 CALL DGEMM('N','N',NCYCLE+1,NCYCLE,NCYCLE+1,ONE,GMTRN,ITERMX+1,
     *           GMRES1,ITERMX+1,ZER,GMRES2,ITERMX+1)
C
      DO ICYCLE= NCYCLE,1,-1
         QTRN= GMTRN(ICYCLE,1)*GAMMA
         DO II= ICYCLE+1,NCYCLE
            QTRN= QTRN-GMRES2(ICYCLE,II)*GMFAC(II)
         ENDDO
         GMFAC(ICYCLE)= QTRN/GMRES2(ICYCLE,ICYCLE)
      ENDDO
C
C     ----- CALCULATE PAI -----
C
      DO MOA= 1,NVIR
         DO MOI= 1,NOCC
            PAI(MOI,MOA)=GMVEC(MOI,MOA,0)
         ENDDO
      ENDDO
      DO ICYCLE= 1,NCYCLE
         DO MOA= 1,NVIR
            DO MOI= 1,NOCC
               PAI(MOI,MOA)= PAI(MOI,MOA)
     *                      +GMFAC(ICYCLE)*GMVEC(MOI,MOA,ICYCLE)
            ENDDO
         ENDDO
      ENDDO
C
 9100 FORMAT(1X,'CYCLE',1X,I2,2X,'Z-VECTOR ERROR=',1X,E13.6)
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGWAI
      SUBROUTINE MPGWAI(NFT16,ENG,PAI,WIJ,WAI,WIJ1,XKIBJ,
     *                  NOCC,NVIR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ENG(*),PAI(NOCC,*),WIJ(NOCC,*),WAI(NOCC,*),WIJ1(NOCC,*),
     *          XKIBJ(NOCC,NOCC,*)
C
      CALL VCLR(WIJ1,1,NOCC*NOCC)
C
C     ----- CALCULATE WAI[II] -----
C
      DO MOA= 1,NVIR
         DO MOI= 1,NOCC
            WAI(MOI,MOA)= WAI(MOI,MOA)-PAI(MOI,MOA)*ENG(MOI)*TWO
         ENDDO
      ENDDO
C
C     ----- CALCULATE A PAI PART OF WIJ[III] -----
C
      NPASS= NVIR/NPROC
      IF(MOD(NVIR,NPROC).NE.0) NPASS= NPASS+1
C
      DO IPASS= 1,NPASS
         IF(GOPARR) THEN
           MOA=(IPASS-1)*NPROC+ME+1
           IF(MOA.GT.NVIR) GO TO 100
         ELSE
           MOA= IPASS
         ENDIF
         ICOUNT=(IPASS-1)*NOCC
         DO MOI= 1,NOCC
            READ(UNIT=NFT16,REC=ICOUNT+MOI)
     *          ((XKIBJ(K,J,MOI),K=1,NOCC),J=1,NOCC)
         ENDDO
         CALL DGEMV('T',NOCC,NOCC*NOCC,-FOUR,XKIBJ,NOCC,
     *              PAI(1,MOA),1,ONE,WIJ1,1)
         DO MOI= 1,NOCC
            PFAC= PAI(MOI,MOA)
            DO J= 1,NOCC
               DO K= 1,NOCC
                  WIJ1(K,J)= WIJ1(K,J)
     *                      +(XKIBJ(K,J,MOI)+XKIBJ(J,K,MOI))*PFAC
               ENDDO
            ENDDO
         ENDDO
      ENDDO
C
  100 IF(GOPARR) CALL DDI_GSUMF(1001,WIJ1,NOCC*NOCC)
      CALL DAXPY(NOCC*NOCC,ONE,WIJ1,1,WIJ,1)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGTAO
      SUBROUTINE MPGTAO(VECO,VECO2,VECV,ENG,PMN,WMN,
     *                  PAI,WIJ,WAB,WAI,WMJ,WMB,
     *                  PAO,WAO,NOCC,NVIR,NBF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION VECO(*),VECO2(NOCC,*),VECV(*),ENG(*),PMN(NBF,*),
     *          WMN(NBF,*),PAI(*),WIJ(*),WAB(*),WAI(*),WMJ(*),WMB(*),
     *          PAO(*),WAO(*)
C
C     ----- TRANSFORM PAI TO PMN -----
C
      CALL DGEMM('N','T',NOCC,NBF,NVIR,TWO,PAI,NOCC,
     *           VECV,NBF,ZER,WMJ,NOCC)
      CALL DGEMM('N','N',NBF,NBF,NOCC,ONE,VECO,NBF,
     *           WMJ,NOCC,-TWO,PMN,NBF)
C
C     ----- TRANSFORM WIJ AND WAI TO WMN -----
C
      CALL DGEMM('N','T',NOCC,NBF,NOCC,ONE,WIJ,NOCC,
     *           VECO,NBF,ZER,WMJ,NOCC)
      CALL DGEMM('N','T',NOCC,NBF,NVIR,ONE,WAI,NOCC,
     *           VECV,NBF,ONE,WMJ,NOCC)
      CALL DGEMM('N','N',NBF,NBF,NOCC,ONE,VECO,NBF,
     *           WMJ,NOCC,ZER,WMN,NBF)
C
C     ----- TRANSFORM WAB TO WMN -----
C
      CALL DGEMM('N','T',NVIR,NBF,NVIR,ONE,WAB,NVIR,
     *           VECV,NBF,ZER,WMB,NVIR)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,VECV,NBF,
     *           WMB,NVIR,ONE,WMN,NBF)
C
      IJ= 0
      DO I= 1,NBF
         DO J= 1,I
            IJ= IJ+1
            PAO(IJ)= PMN(I,J)+ PMN(J,I)
            WAO(IJ)= WMN(I,J)+ WMN(J,I)
         ENDDO
      ENDDO
C
      NBF3=(NBF*NBF+NBF)/2
      CALL DSCAL(NBF3,0.5D+0,PAO,1)
      CALL DSCAL(NBF3,0.5D+0,WAO,1)
C
C     ----- ADD HF+MP2 ENERGY-WEIGHTED DENSITY -----
C
      IJ= 0
      DO I= 1,NBF
         DO J= 1,I
            IJ= IJ+1
            DUM= ZER
            DO K= 1,NOCC
               DUM= DUM-ENG(K)*VECO2(K,I)*VECO2(K,J)
            ENDDO
            WAO(IJ)= WAO(IJ)+DUM*TWO
         ENDDO
      ENDDO
      CALL DAWRIT(IDAF,IODA,WAO,NBF3,309,0)
C
C     ----- SAVE PMN ON 307 -----
C
      CALL DAWRIT(IDAF,IODA,PAO,NBF3,307,0)
C
C     ----- SAVE PHF ON 308 -----
C
      CALL DAREAD(IDAF,IODA,WAO,NBF3,16,0)
      CALL DAWRIT(IDAF,IODA,WAO,NBF3,308,0)
C
C     ----- FORM AND SAVE TOTAL HF+MP2 DENSITY -----
C
      CALL DAXPY(NBF3,ONE,WAO,1,PAO,1)
      CALL DAWRIT(IDAF,IODA,PAO,NBF3,16,0)
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGJKD
      SUBROUTINE MPGJKD(NFT20,VECO,XINT,TIJ,TMJ,TMN,PHF,PMP2,
     *                  DE1,NLOOP,ILOOP,NOCC,NBF,NOAC,NANGM,NEEDGRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL POPLE,HONDO,MP2
      LOGICAL PACK2E,OUT,DBG,GOPARR,DSKWRK,MASWRK
C
      INTEGER IPERM(4), INVPERM(4)
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000, MXAO=8192)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (TEN=10.0D+00, ONE=1.0D+00, ZER=0.0D+00)
      PARAMETER (TENM9=1.0D-09, TENM11=1.0D-11)
      PARAMETER (TENM20=1.0D-20, HALF=0.5D+00, TENM12=1.0D-12)
C
      COMMON /FMCOM / X(1)
      COMMON /DSHLNO/ LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJD,KLD,IJ,KL
      COMMON /DSHLT / RTOL,DTOL,VTOL1,VTOL2,VTOLS,OUT,DBG
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
C     COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SHLBAS/ MAXTYP,MAXNUM
      COMMON/DERMEM/IWFN,IXCH,INIJG,IGINT,IFINT,ISINT,IIJKLG,
     1 IDAB,ICHRG,IXY,IXZ,IYZ,IX,IY,IZ,ISJ,ISK,ISL,IGIJKL,IGNKL,IGNM,
     2 IDIJ,IDKL,IB00,IB01,IB10,IC00,ID00,IF00,
     3 IAAI,IAAJ,IBBK,IBBL,IFI,IFJ,IFK,IFL,
     4 ISII,ISJJ,ISKK,ISLL,ISIJ,ISIK,ISIL,ISJK,ISJL,ISKL,
     5 IDIJSI,IDIJSJ,IDKLSK,IDKLSL,IABV,ICV,IRW
      COMMON/INDD80/IMAX,JMAX,KKKMAX,LMAX
C
      DIMENSION VECO(*),XINT(*),TIJ(*),TMJ(*),TMN(*),PHF(*),PMP2(*),
     *          DE1(*),NLOOP(0:NPROC-1,4),ILOOP(NSHELL,0:NPROC-1,2)
C
      DATA CHECK/8HCHECK   /
C
C       THE DRIVER FOR THE TWO ELECTRON GRADIENT
C
C              CALCULATE THE LARGEST SHELL TYPE
C
      NBF3=(NBF*NBF+NBF)/2
      CALL BASCHK(MAXTYP)
      IF(MAXTYP.GT.5) THEN
         IF(MASWRK) WRITE(IW,*)
     *        'MP2 GRADIENT INTEGRALS PROGRAMMED FOR SPDFG ONLY'
         CALL ABRT
      END IF
C
C     ----- SET STARTING PARAMETERS -----
C
      MP2   = .TRUE.
      HONDO = .TRUE.
      POPLE = .TRUE.
      IF(IGRDTYP.EQ.2)          POPLE = .FALSE.
      IF(POPLE.AND.MAXTYP.LT.2) HONDO = .FALSE.
C
C      CUTOFF IS THE SCHWARZ SCREENING CUT OFF
C      DABCUT IS THE TWO PARTICLE DENSITY CUT OFF
C
      CUTOFF=TENM9
      IF(.NOT.POPLE) CUTOFF=CUTOFF/TEN
      CUTOFF2=CUTOFF*HALF
C
      ZBIG = 0.0D+00
      DO ISH=1,NSHELL
         I1=KSTART(ISH)
         I2=I1+KNG(ISH)-1
         DO IG=I1,I2
            IF(EX(IG).GT.ZBIG) ZBIG = EX(IG)
         ENDDO
      ENDDO
      DABCUT=TENM11
      IF(ZBIG.GT.1.0D+06) DABCUT = DABCUT/TEN
      IF(ZBIG.GT.1.0D+07) DABCUT = DABCUT/TEN
C
C      VTOLS ARE CUT OFFS USED BY THE POPLE PACKAGE
C      CURRENT VALUES ARE FROM HONDO 8, SEE G92 FOR OTHER POSSIBILITIES
C
      Q4    = ONE
      VTOL1= TENM12
      VTOL2= TENM12
      VTOLS= TENM20
      DTOL = TEN**(-ITOL)
      RTOL = RLN10*ITOL
C
C      INITIALIZE THE INTEGRAL BLOCK COUNTERS TO ZERO
C
      IISKIP1= 0
      IISKIP2= 0
      IDID   = 0
      LENGTH = 1
      NC     = 1
      LDF    = 1
      ICOUNT = 0
C
C  SET UP MEMORY POINTERS FOR DERIVATIVE INTEGRALS
C
      CALL VALFMP(LOADFM)
      CALL JKDMEM(1,LOADFM,IADDR,LENGTH,
     *            MINXYZ,MAXXYZ,MINVEC,POPLE,MP2)
      NEED = IADDR - LOADFM
      CALL GETFMP(NEED)
C
      CALL VALFMP(LOADFM)
      IDAB = 1    + LOADFM
      LAST = IDAB + MAXTYP**4
      NEED2= LAST - LOADFM -1
      CALL GETFMP(NEED2)
      NEEDGRD = NEED + NEED2
C
C  SAVE 1-PARTICLE GRADIENT IN TEMPORARY BUFFER, ZERO 2E- GRADIENT
C
      CALL DCOPY(NAT*3,DE,1,DE1,1)
      CALL DCOPY(NAT*3,ZER,0,DE,1)
C
      IF(EXETYP.EQ.CHECK) GO TO 200
C
      CALL DAREAD(IDAF,IODA,PHF,NBF3,308,0)
      CALL DAREAD(IDAF,IODA,PMP2,NBF3,307,0)
C
C  PRE-COMPUTE FACTORS FOR INTEGRAL EVALUATION
C
      IF(HONDO) CALL OEDHND(X(INIJG),X(ICHRG))
      DTOL = DTOL*DTOL
C
C  I-SHELL
C
      IEND= NSHELL
      IF(GOPARR) IEND= NLOOP(ME,3)
      DO I1= 1,IEND
        II= I1
        IF(GOPARR) II= ILOOP(I1,ME,2)
        ISHL = KMAX(II)-KMIN(II)+1
C
C  K-SHELL
C
        DO KK= 1,II
          KSHL = KMAX(KK)-KMIN(KK)+1
          NMINK= KLOC(KK)-1
C
C  THIRD BACK-TRANSFORMATION OF AMPLITUDES
C
          CALL MPGBT3(NFT20,VECO,TIJ,TMJ,NOCC,NBF,NOAC,ISHL,ICOUNT,KSHL,
     *                NMINK,NANGM)
C
C  J-SHELL
C
          DO JJ= 1,NSHELL
            JSHL = KMAX(JJ)-KMIN(JJ)+1
            NMINJ= KLOC(JJ)-1
C
C  L-SHELL
C
            DO LL= 1,NSHELL
              LSHL= KMAX(LL)-KMIN(LL)+1
              NMINL= KLOC(LL)-1
              ISH = II
              JSH = JJ
              KSH = KK
              LSH = LL
C
C  SINGLE-CENTRE NUCLEAR GRADIENTS ARE NULL
C
              IF(KATOM(JSH).EQ.KATOM(ISH).AND.
     *           KATOM(KSH).EQ.KATOM(ISH).AND.
     *           KATOM(LSH).EQ.KATOM(ISH)) GO TO 100
C
C  SCHWARZ INEQUALITY TEST
C
              IIJJ= IA(MAX0(II,JJ))+MIN0(II,JJ)
              KKLL= IA(MAX0(KK,LL))+MIN0(KK,LL)
              GMAX= XINT(IIJJ)*XINT(KKLL)
              IF(GMAX.LT.CUTOFF) THEN
                IISKIP1= IISKIP1+1
                GO TO 100
              ENDIF
C
C  USE POPLE CODE FOR ANY PURE SP INTEGRAL BLOCKS,
C  USE HONDO RYS POLYNOMIAL CODE FOR OTHER BLOCKS
C
              POPLE = .TRUE.
              IF(IGRDTYP.EQ.2)    POPLE=.FALSE.
              IF(KTYPE(ISH).GT.2) POPLE=.FALSE.
              IF(KTYPE(JSH).GT.2) POPLE=.FALSE.
              IF(KTYPE(KSH).GT.2) POPLE=.FALSE.
              IF(KTYPE(LSH).GT.2) POPLE=.FALSE.
C
C  ENSURE SHELL LABELS ARE IN CANONICAL ORDER
C
              INEW = ISH
              JNEW = JSH
              KNEW = KSH
              LNEW = LSH
              IPERM(1)= 1
              IPERM(2)= 2
              IPERM(3)= 3
              IPERM(4)= 4
              IF(INEW.LT.JNEW) THEN
                ITMP= INEW
                INEW= JNEW
                JNEW= ITMP
                ITMP= IPERM(1)
                IPERM(1)= IPERM(2)
                IPERM(2)= ITMP
              ENDIF
              IF(KNEW.LT.LNEW) THEN
                ITMP= KNEW
                KNEW= LNEW
                LNEW= ITMP
                ITMP= IPERM(3)
                IPERM(3)= IPERM(4)
                IPERM(4)= ITMP
              ENDIF
              IF(IIJJ.LT.KKLL) THEN
                ITMP= INEW
                INEW= KNEW
                KNEW= ITMP
                ITMP= IPERM(1)
                IPERM(1)= IPERM(3)
                IPERM(3)= ITMP
                ITMP= JNEW
                JNEW= LNEW
                LNEW= ITMP
                ITMP= IPERM(2)
                IPERM(2)= IPERM(4)
                IPERM(4)= ITMP
                ITMP= IIJJ
                IIJJ= KKLL
                KKLL= ITMP
              ENDIF
C
C      ----- POPLE SET UP -----
C
              IF(POPLE) THEN
C
C  NON-CANONICAL 4-LABEL SYMMETRY WEIGHT
C
                SFAC = 0.125D+00
                IF (KSH.NE.ISH) SFAC = 0.25D+00
C
C  ENSURE SHELL ANGULAR MOMENTA ARE IN CANONICAL ORDER
C
                IMAX=KTYPE(INEW)-1
                JMAX=KTYPE(JNEW)-1
                KKKMAX=KTYPE(KNEW)-1
                LMAX=KTYPE(LNEW)-1
                IF(IMAX.LT.JMAX) THEN
                  ITMP = INEW
                  INEW = JNEW
                  JNEW = ITMP
                  ITMP = IPERM(1)
                  IPERM(1) = IPERM(2)
                  IPERM(2) = ITMP
                END IF
                IF(KKKMAX.LT.LMAX) THEN
                  ITMP = KNEW
                  KNEW = LNEW
                  LNEW = ITMP
                  ITMP = IPERM(3)
                  IPERM(3) = IPERM(4)
                  IPERM(4) = ITMP
                END IF
                IF((IMAX+JMAX).LT.(KKKMAX+LMAX)) THEN
                  ITMP = INEW
                  INEW = KNEW
                  KNEW = ITMP
                  ITMP = IPERM(1)
                  IPERM(1) = IPERM(3)
                  IPERM(3) = ITMP
                  ITMP = JNEW
                  JNEW = LNEW
                  LNEW = ITMP
                  ITMP = IPERM(2)
                  IPERM(2) = IPERM(4)
                  IPERM(4) = ITMP
                  ITMP = IIJJ
                  IIJJ = KKLL
                  KKLL = ITMP
                END IF
C
C   GENERATE THE INVERSE PERMUTATION OF IPERM
C   TO MAP THE NON-SEPARABLE INDICES ONTO THE
C   PERMUTED SHELL ORDERING
C
                DO I= 1,4
                  DO J= 1,4
                    IF(IPERM(J).EQ.I) INVPERM(I)= J
                  END DO
                END DO
                IMAX=3*(KTYPE(INEW)-1)+1
                JMAX=3*(KTYPE(JNEW)-1)+1
                KKKMAX=3*(KTYPE(KNEW)-1)+1
                LMAX=3*(KTYPE(LNEW)-1)+1
                JTYPE=(IMAX+JMAX+KKKMAX+KKKMAX+LMAX-2)/3
                IAT = KATOM(INEW)
                JAT = KATOM(JNEW)
                KAT = KATOM(KNEW)
                LAT = KATOM(LNEW)
              ELSE
C
C    ----- HONDO SET UP -----
C
C  SHELL 4-LABEL WEIGHT IN THE NON-CANONICAL LIST
C
                SFAC = ONE
                IF(IIJJ.NE.KKLL) SFAC = HALF
                IF(ISH.NE.JSH) SFAC = SFAC*HALF
                IF(KSH.NE.LSH) SFAC = SFAC*HALF
                IF(KSH.NE.ISH) SFAC = SFAC*2.0D+00
C
C   GENERATE INVERSE PERMUTATION
C
                DO I = 1, 4
                  DO J = 1, 4
                    IF (IPERM(J).EQ.I) INVPERM(I) = J
                  END DO
                END DO
C
C  SET POINTERS TO THE IJ AND KL CHARGE DISTRIBUTIONS
C
                CALL OEDRD(X(INIJG),NIJ,NIJ0,IIJJ)
                CALL OEDRD(X(INIJG),NKL,NKL0,KKLL)
C
C  SELECT CENTERS FOR DERIVATIVES
C
                CALL JKDATM(INEW,JNEW,KNEW,LNEW)
C
C  SET INDICES FOR SHELL BLOCK
C
                CALL JKDSHL(INEW,JNEW,KNEW,LNEW)
                CALL JKDNDX(X(IIJKLG))
              END IF
C
C  FOURTH BACK-TRANSFORMATION OF AMPLITUDES
C
              CALL MPGBT4(VECO,TMJ,TMN,ISHL,JSHL,KSHL,LSHL,NMINJ,NMINL,
     *                    NOCC,NBF,NOAC,NANGM)
C
C  COMBINE TERMS OF THE 2-PARTICLE DENSITY MATRIX
C
              CALL MPG2PDM(TMN,PHF,PMP2,X(IDAB),INEW,JNEW,KNEW,LNEW,
     *                     INVPERM,SFAC,DABMAX,POPLE,NANGM)
C
C   FINE SCREENING, ON INTEGRAL VALUE TIMES DENSITY FACTOR
C
              IF(DABMAX*GMAX.LT.CUTOFF2) THEN
                IISKIP2 = IISKIP2+1
                GO TO 100
              END IF
C
C  CONTRACT DENSITY WITH THE DERIVATIVE INTEGRALS
C  AND ADD DIRECTLY TO GRADIENT VECTOR
C
              IDID = IDID+1
              IF (POPLE) THEN
                CALL JKDG80(DABMAX,INEW,JNEW,KNEW,LNEW,
     *                      JTYPE,IAT,JAT,KAT,LAT)
              ELSE
            CALL JKDSPD(NIJ0,NKL,NKL0,X(ICHRG),
     *      X(IGINT),X(IFINT),X(ISINT),X(IIJKLG),X(IGIJKL),
     1      X(IGNKL),X(IGNM),X(IXY),X(IXZ),X(IYZ),X(IX),X(IY),X(IZ),
     2      X(ISJ),X(ISK),X(ISL),X(IB00),X(IB01),X(IB10),X(IC00),
     3      X(ID00),X(IF00),X(IDIJ),X(IDKL),X(IDIJSI),X(IDIJSJ),
     4      X(IDKLSK),X(IDKLSL),X(IABV),X(ICV),X(IRW),X(IAAI),X(IAAJ),
     5      X(IBBK),X(IBBL),X(IFI),X(IFJ),X(IFK),X(IFL),X(ISII),X(ISJJ),
     6      X(ISKK),X(ISLL),X(ISIJ),X(ISIK),X(ISIL),X(ISJK),X(ISJL),
     7      X(ISKL),X(IDAB),MAXXYZ,FC,NC,DF,LDF,NBF,DDA,Q4,MINVEC,
     8      DABCUT,DABMAX)
              END IF
C
  100         CONTINUE
            ENDDO
          ENDDO
        ENDDO
        ICOUNT= ICOUNT+ISHL
      ENDDO
C
  200 CONTINUE
      CALL RETFMP(NEED2)
      CALL RETFMP(NEED)
C
      IF(GOPARR) THEN
        CALL DDI_GSUMI(1909,IISKIP1,1)
        CALL DDI_GSUMI(1908,IISKIP2,1)
        CALL DDI_GSUMI(1907,IDID   ,1)
        CALL DDI_GSUMF(1109,DE,NAT*3)
      ENDIF
      IF(MASWRK  .AND.  EXETYP.NE.CHECK)
     *      WRITE(IW,9999) IISKIP1,IISKIP2,IDID
C
      CALL DAXPY(NAT*3,ONE,DE1,1,DE,1)
C
      CALL DFINAL(1)
      IREST = 0
      RETURN
C
 9999 FORMAT(1X,'THE COARSE/FINE SCHWARZ SCREENINGS SKIPPED ',I10,'/',
     *          I10,' BLOCKS.'/
     *     1X,'THE NUMBER OF GRADIENT INTEGRAL BLOCKS COMPUTED WAS',I10)
      END
C
C*MODULE MP2GR2  *DECK MPG2PDM
      SUBROUTINE MPG2PDM(TMN,PHF,PMP2,DAB,INEW,JNEW,KNEW,LNEW,
     *                   INVPERM,SFAC,DABMAX,POPLE,NANGM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL POPLE
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (ZER=0.0D+00, PT5=0.5D+00, F04=4.0D+00)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CCG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLLMN/ IGXYZ(4,35),JGXYZ(4,35),KGXYZ(4,35),LGXYZ(4,35)
      COMMON /SHLNRM/ PNRM(84)
      COMMON /GSPG80/ E(4,4,4,4)
C
      DIMENSION TMN(NANGM,NANGM,NANGM,*),PHF(*),PMP2(*),DAB(*)
      DIMENSION INVPERM(4),IINDX(4)
C
C  COMBINE TERMS OF 2-PARTICLE DENSITY MATRIX FOR A SHELL OF
C  DERIVATIVE INTEGRALS - COMPUTING SCF AND SEPARABLE TERMS
C
      DABMAX= ZER
      MINI= KMIN(INEW)
      MINJ= KMIN(JNEW)
      MINK= KMIN(KNEW)
      MINL= KMIN(LNEW)
      MAXI= KMAX(INEW)
      MAXJ= KMAX(JNEW)
      MAXK= KMAX(KNEW)
      MAXL= KMAX(LNEW)
      LOCI= KLOC(INEW)-MINI
      LOCJ= KLOC(JNEW)-MINJ
      LOCK= KLOC(KNEW)-MINK
      LOCL= KLOC(LNEW)-MINL
C
      IF(POPLE) THEN
         DO I=1,4
            IGXYZ(1,I)=64*(I-1)+1
            JGXYZ(1,I)=16*(I-1)
            KGXYZ(1,I)= 4*(I-1)
            LGXYZ(1,I)=   (I-1)
         ENDDO
C
         DO L=MINL,MAXL
            IINDX(4) = L-(MINL-1)
            NNU = LOCL+L
            DO K=MINK,MAXK
               IINDX(3) = K-(MINK-1)
               NMU = LOCK+K
               MUNU= IA(MAX0(NMU,NNU))+MIN0(NMU,NNU)
               DO J=MINJ,MAXJ
                  IINDX(2) = J-(MINJ-1)
                  NSI = LOCJ+J
                  DO I=MINI,MAXI
                     IINDX(1) = I-(MINI-1)
                     NLA = LOCI+I
                     LASI= IA(MAX0(NLA,NSI))+MIN0(NLA,NSI)
                     MUSI= IA(MAX0(NMU,NSI))+MIN0(NMU,NSI)
                     LANU= IA(MAX0(NLA,NNU))+MIN0(NLA,NNU)
                     MULA= IA(MAX0(NMU,NLA))+MIN0(NMU,NLA)
                     NUSI= IA(MAX0(NNU,NSI))+MIN0(NNU,NSI)
C
                     DF1=(PHF(MUNU)+PMP2(MUNU))*PHF(LASI)+
     *                    PHF(MUNU)*PMP2(LASI)
                     DQ1=(PHF(MUSI)+PMP2(MUSI))*PHF(LANU)+
     *                    PHF(MUSI)*PMP2(LANU)+
     *                   (PHF(MULA)+PMP2(MULA))*PHF(NUSI)+
     *                    PHF(MULA)*PMP2(NUSI)
C
C  MP2 NON-SEPARABLE TERM
C
                     IP = IINDX(INVPERM(1))
                     JP = IINDX(INVPERM(2))
                     KP = IINDX(INVPERM(3))
                     LP = IINDX(INVPERM(4))
                     DF1=-DQ1+F04*(DF1+F04*TMN(LP,JP,KP,IP))
                     DF1= DF1*SFAC
                     IF(DABMAX.LT.ABS(DF1)) DABMAX= ABS(DF1)
                     E(I,J,K,L)= DF1
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ELSE
C
C D AND HIGHER FUNCTIONS OR HONDO ONLY RUNS
C
         DO I=MINI,MAXI
            P1I= PNRM(I)
            IINDX(1) = I-(MINI-1)
            NLA = LOCI+I
            DO J=MINJ,MAXJ
               P2J= P1I*PNRM(J)
               IINDX(2) = J-(MINJ-1)
               NSI = LOCJ+J
               LASI= IA(MAX0(NLA,NSI))+MIN0(NLA,NSI)
               DO K=MINK,MAXK
                  P3K= P2J*PNRM(K)
                  IINDX(3) = K-(MINK-1)
                  NMU = LOCK+K
                  DO L=MINL,MAXL
                     P4L= P3K*PNRM(L)
                     IINDX(4) = L-(MINL-1)
                     NNU = LOCL+L
                     MUNU= IA(MAX0(NMU,NNU))+MIN0(NMU,NNU)
                     MUSI= IA(MAX0(NMU,NSI))+MIN0(NMU,NSI)
                     LANU= IA(MAX0(NLA,NNU))+MIN0(NLA,NNU)
                     MULA= IA(MAX0(NMU,NLA))+MIN0(NMU,NLA)
                     NUSI= IA(MAX0(NNU,NSI))+MIN0(NNU,NSI)
C
                     DF1=(PHF(MUNU)+PMP2(MUNU))*PHF(LASI)+
     *                    PHF(MUNU)*PMP2(LASI)
                     DQ1=(PHF(MUSI)+PMP2(MUSI))*PHF(LANU)+
     *                    PHF(MUSI)*PMP2(LANU)+
     *                   (PHF(MULA)+PMP2(MULA))*PHF(NUSI)+
     *                    PHF(MULA)*PMP2(NUSI)
C
C  MP2 NON-SEPARABLE TERM
C
                     IP = IINDX(INVPERM(1))
                     JP = IINDX(INVPERM(2))
                     KP = IINDX(INVPERM(3))
                     LP = IINDX(INVPERM(4))
                     DF1=-DQ1+F04*(DF1+F04*TMN(LP,JP,KP,IP))
                     DF1= DF1*P4L*SFAC
                     IF(NMU .EQ.NNU ) DF1= DF1*PT5
                     IF(NLA .EQ.NSI ) DF1= DF1*PT5
                     IF(MUNU.EQ.LASI) DF1= DF1*PT5
                     IF(DABMAX.LT. ABS(DF1)) DABMAX = ABS(DF1)
                     IJKL=IGXYZ(1,I)+JGXYZ(1,J)+KGXYZ(1,K)+LGXYZ(1,L)
                     DAB(IJKL)= DF1
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
C
      IF(INEW.EQ.JNEW) THEN
         DO L=MINL,MAXL
            DO K=MINK,MAXK
               DO J=MINJ,MAXJ
                  DO I=MINI,MAXI
                     IF(POPLE) THEN
                        DF1=(E(I,J,K,L)+E(J,I,K,L))*PT5
                        E(I,J,K,L)= DF1
                        E(J,I,K,L)= DF1
                     ELSE
                        IJKL=IGXYZ(1,I)+JGXYZ(1,J)+KGXYZ(1,K)+LGXYZ(1,L)
                        JIKL=IGXYZ(1,J)+JGXYZ(1,I)+KGXYZ(1,K)+LGXYZ(1,L)
                        DF1=(DAB(IJKL) + DAB(JIKL))*PT5
                        DAB(IJKL) = DF1
                        DAB(JIKL) = DF1
                     END IF
                  END DO
               END DO
            END DO
         END DO
      END IF
C
      IF(KNEW.EQ.LNEW) THEN
         DO L=MINL,MAXL
            DO K=MINK,MAXK
               DO J=MINJ,MAXJ
                  DO I=MINI,MAXI
                     IF(POPLE) THEN
                        DF1=(E(I,J,K,L)+E(I,J,L,K))*PT5
                        E(I,J,K,L)= DF1
                        E(I,J,L,K)= DF1
                     ELSE
                        IJKL=IGXYZ(1,I)+JGXYZ(1,J)+KGXYZ(1,K)+LGXYZ(1,L)
                        IJLK=IGXYZ(1,I)+JGXYZ(1,J)+KGXYZ(1,L)+LGXYZ(1,K)
                        DF1=(DAB(IJKL) + DAB(IJLK))*PT5
                        DAB(IJKL) = DF1
                        DAB(IJLK) = DF1
                     END IF
                  END DO
               END DO
            END DO
         END DO
      END IF
C
      IJ = IA(MAX0(INEW,JNEW))+MIN0(INEW,JNEW)
      KL = IA(MAX0(KNEW,LNEW))+MIN0(KNEW,LNEW)
      IF(IJ.EQ.KL) THEN
         DO L=MINL,MAXL
            DO K=MINK,MAXK
               DO J=MINJ,MAXJ
                  DO I=MINI,MAXI
                     IF(POPLE) THEN
                        DF1=(E(I,J,K,L)+E(K,L,I,J))*PT5
                        E(I,J,K,L)= DF1
                        E(K,L,I,J)= DF1
                     ELSE
                        IJKL=IGXYZ(1,I)+JGXYZ(1,J)+KGXYZ(1,K)+LGXYZ(1,L)
                        KLIJ=IGXYZ(1,K)+JGXYZ(1,L)+KGXYZ(1,I)+LGXYZ(1,J)
                        DF1=(DAB(IJKL) + DAB(KLIJ))*PT5
                        DAB(IJKL) = DF1
                        DAB(KLIJ) = DF1
                     END IF
                  END DO
               END DO
            END DO
         END DO
      END IF
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGBT3
      SUBROUTINE MPGBT3(NFT20,VECO,TIJ,TMJ,NOCC,NBF,NOAC,ISHL,ICOUNT,
     *                  KSHL,NMINK,NANGM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00, ZER=0.0D+00)
C
      DIMENSION VECO(*),TIJ(NOAC*NOAC,NANGM,*),TMJ(NOAC*NBF,NANGM,*)
C
C     3RD-QUARTER BACK-TRANSFORMATION OF AMPLITUDES
C
      NOAC2= NOAC*NOAC
      DO II= 1,ISHL
         IREC=(ICOUNT+II-1)*NBF+NMINK
         DO KK= 1,KSHL
            READ(UNIT=NFT20,REC=IREC+KK) (TIJ(IJ,KK,II),IJ=1,NOAC2)
         ENDDO
      ENDDO
C
      DO II= 1,ISHL
         DO KK= 1,KSHL
            CALL DGEMM('N','N',NOAC,NBF,NOAC,ONE,TIJ(1,KK,II),NOAC,
     *                 VECO,NOCC,ZER,TMJ(1,KK,II),NOAC)
         ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE MP2GR2  *DECK MPGBT4
      SUBROUTINE MPGBT4(VECO,TMJ,TMN,ISHL,JSHL,KSHL,LSHL,NMINJ,NMINL,
     *                  NOCC,NBF,NOAC,NANGM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VECO(NOCC,*),TMJ(NOAC,NBF,NANGM,*),
     *          TMN(NANGM,NANGM,NANGM,*)
C
C     4TH-QUARTER BACK-TRANSFORMATION OF AMPLITUDES
C
      DO I= 1,ISHL
         DO K= 1,KSHL
            DO J= 1,JSHL
               JJ= NMINJ+J
               DO L= 1,LSHL
                  LL= NMINL+L
                  TMN(L,J,K,I)= DDOT(NOAC,VECO(1,LL),1,TMJ(1,JJ,K,I),1)
               ENDDO
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
