C 21 APR 08 - MWS - INCLUDE EXETYP=CHECK OPTION
C  4 MAR 08 - KI  - ALLOW ABELIAN GROUP SYMMETRY
C 20 AUG 07 - DGF - SYNCHRONISE FMORUN
C 14 FEB 07 - MWS - USE A WRAPPER NAME FOR GET ENVIRONMENT CALL
C 27 DEC 06 - KI  - PARALLEL PROGRAM FROM IMS
C
C*MODULE MP2IMS  *DECK MP2IMS
      SUBROUTINE MP2IMS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXSH=5000, MXGTOT=20000)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /FMCOM / XX(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
      CALL VALFM(LOADFM)
      LNUMP  = LOADFM + 1
      LNLOOP = LNUMP  + NPROC + 1
      LILOOP = LNLOOP + NPROC*2
      LIJEND = LILOOP + NSHELL*NPROC
      LILEN  = LIJEND + NPROC
      LIRECV = LILEN  + NPROC
      LAST   = LIRECV + NPROC
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      NCPU = NPROC
      CALL MP2IMSX(XX(LNUMP),XX(LNLOOP),XX(LILOOP),XX(LIJEND),
     *             XX(LILEN),XX(LIRECV),NCPU)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE MP2IMS  *DECK MP2IMSX
      SUBROUTINE MP2IMSX(NUMP,NLOOP,ILOOP,IJEND,ILEN,IRECV,NCPU)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ISGDDI,PAROUT,INITGDDI,ABEL,ABELPT,
     *        GOPARR,DSKWRK,MASWRK,DSKSAV,DLB,SLB,DIVANG
      PARAMETER (MXUNIT=299)
      CHARACTER*1   NULL
      CHARACTER*256 PATHNM,ENVBUF
      COMMON /ENVIR/ ENVBUF(MXUNIT)
      COMMON /OPNNFT/ NFTOPN(MXUNIT),NODEXT(MXUNIT),IOSMP(2)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXRT=100, MXATM=2000)
      PARAMETER (MXAO=8192)
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, P12=1.2D+00)
C
C        NOTE THAT -NCPU- IS AN ALIAS FOR -NPROC-
      DIMENSION NUMP(0:NCPU),NLOOP(2,NCPU),ILOOP(NSHELL,NCPU),
     *          IJEND(NCPU),ILEN(NCPU),IRECV(NCPU)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /DIRMEM/ LPQRJ,IVEC,IDEN,IFCK,NUM,NOCC,NDIM,MINPQ,MAXPQ
      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 /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,NPGLOB,NNGLOB,JBTYP,
     *                ISGDDI,PAROUT,INITGDDI
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2DAF/ IDAF20,NAV20
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      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 /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
      DATA CHECK/8HCHECK   /,ANONE/8HNONE    /
C
C     ----- MAIN CLOSED-SHELL SEMI-DIRECT MP2 ENERGY DRIVER -----
C
      DSKSAV= DSKWRK
      DSKWRK= .TRUE.
      DLB = IBTYP.EQ.1
      SLB = .NOT.DLB
C
C        THE INTEGRAL TRANSFORMATION CAN EXPLOIT ABELIAN GROUPS
C
      ABEL=ABELPT()
      IF(.NOT.ABEL) CALL SYMOFF
C
      E2S= ZER
      E2T= ZER
      CUTOFF=TOL
      QQ4=ONE
C
C     SUBTRACT FMO PROJECTED OUT ORBITALS.
C
      IF(NFG.NE.0) NO = NO - NORBPROJ
C
      IF(MASWRK) WRITE(IW,9000) NOA,NOA-NACORE,NO-NOA,NBF,NSHELL,CUTOFF
      IF(MASWRK.AND.NFG.NE.0.AND.NORBPROJ.NE.0) WRITE(IW,9010) NORBPROJ
C
      CALL SETFMPI(MEMPRI)
C
C    ----- VARIABLES FOR MP2 ENERGY -----
C
C NO       : NUMBER OF CANONICAL MOLECULAR ORBITALS
C NBF      : NUMBER OF BASIS FUNCTIONS
C NSHELL   : NUMBER OF SHELLS
C NOA      : NUMBER OF OCCUPIED ORBITALS
C NACORE   : NUMBER OF CORE ORBITALS
C NOC      : NUMBER OF CORRELATED ORBITALS
C NVIR     : NUMBER OF VIRTUAL ORBITALS
C NORB     : NUMBER OF ACTIVE ORBITALS
C NOC2     : NUMBER OF IJ PAIR
C NIJ3     : NUMBER OF IJ PAIR IN A BATCH AT 3RD TRANS.
C NIJ4     : NUMBER OF IJ PAIR IN A BATCH AT 4TH TRANS.
C
      NOC  = NOA - NACORE
      NVIR = NO  - NOA
      NOC2 =(NOC*NOC+NOC)/2
      NBF2 = NBF*NBF
      NBF3 =(NBF*NBF+NBF)/2
      NSH2 =(NSHELL*NSHELL+NSHELL)/2
C
      NUM  = NBF
      IDEN = 0
      IFCK = 0
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    ----- SIZE OF DISK CACHE -----
C
      NCACHE= 1000000
      NIJ3 = NCACHE/NVIR
      IF(NIJ3.GT.NOC2) NIJ3= NOC2
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    ----- MEMORY POINTERS FOR MP2 ENERGY -----
C
C IVEC    : MATRIX OF MO COEFFICIENTS        NBF*NBF
C IVECO   : TRANSPOSED MATRIX OF             NOC*NBF
C           OCCUPIED MO COEEFICIENTS
C IVECV   : TRANSPOSED MATRIX OF             NVIR*NBF
C           VIRTUAL MO COEEFICIENTS
C IENG    : ORBITAL ENERGIES                 NBF
C ILAB    : MOSLAB    |                      NBF
C IIRP    : MOSIRP    |---> SEE SYM          NBF
C IDEG    : MOSDEG    |                      NBF
C IGHND   : AO INTEGRALS FOR D,F,G           MAXG --> NANGM**4
C IXINT   : EXCHANGE INTEGRALS               (NSHELL**2 + NSHELL)/2
C
C IINT    : AO INTEGRALS                     (NANGM**3)*NBF
C IPTR1   : PARTIALLY TRANSFORMED INTEGRAL   (NANGM**3)*NOC
C IPTR2   : PARTIALLY TRANSFORMED INTEGRAL   NANGM*NOC2*NBF
C IPTR3   : PARTIALLY TRANSFORMED INTEGRAL   NANGM*NVIR*NIJ3
C
C IDDIS   : WORK SPACE FOR SEND              
C                                  NIJ4*NVIR*MXP*NPROC(~=NIJ4*NVIR*NBF)
C JPTR3   : PARTIALLY TRANSFORMED INTEGRAL   NIJ4*NVIR*NBF
C IENGO   : I+J MO ENERGIES                  NOC2
C
      CALL VALFMP(LOADFM)
      IVEC  = 1     + LOADFM
      IVECO = IVEC  + NBF2
      IVECV = IVECO + NOC*NBF
      IVECV2= IVECV + NVIR*NBF
      IENG  = IVECV2+ NVIR*NBF
      ILAB  = IENG  + NBF
      IIRP  = ILAB  + NBF
      IDEG  = IIRP  + NBF
      INSYM = IDEG  + NBF
      IISYM = INSYM + NT
      LAST  = IISYM + NT*NVIR
      NEEDA = LAST  - LOADFM -1
      CALL GETFMP(NEEDA)
C
      CALL VALFMP(LOADFM)
      IGHND = 1     + LOADFM
      IXINT = IGHND + MAXG
      IINT  = IXINT + NSH2
      IPTR1 = IINT  +(NANGM**3)*NBF
      IPTR3 = IPTR1 +(NANGM**3)*NOC
      IWORK = IPTR3 + NVIR*NIJ3
      LAST  = IWORK + NBF*NIJ3
      NEEDB = LAST  - LOADFM -1
      CALL GETFMP(NEEDB)
C
C     ----- DISTRIBUTE P OF (PQ|RS) -----
C     ----- STATIC LOAD BALANCING -----
C
      MXP=NBF/NPROC
      IF(SLB .AND. GOPARR) THEN
         IPROC= 1
         IFUNC= NBF/NPROC
         NUMP(0)= 0
         DO I=1,NPROC
            NUMP(I)= NSHELL
            NLOOP(2,I)= 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(2,IPROC)=MXP
         ENDIF
         DO I=1,IPROC-1
            IF(NUMP(I).NE.NSHELL) THEN
               IMXP= KLOC(NUMP(I)+1)-KLOC(NUMP(I-1)+1)
               NLOOP(2,I)=IMXP
               IF(IMXP.GT.MXP) MXP= IMXP
            ELSE
               NLOOP(2,I)=0
            ENDIF
         ENDDO
C
C            DYNAMIC LOAD BALANCING SETS REAL NLOOP VALUES IN -MP2SDI-
C            SERIAL DOESN'T SEEM TO NEED TO HAVE NLOOP SET TO ANYTHING
C
      ELSE
         DO IPROC=1,NPROC
            DO I=1,2
               NLOOP(I,IPROC) = 0
            ENDDO
         ENDDO
      ENDIF
C
C     ----- PRINT OUT NEEDED MEMORY AND DISK SPACE -----
C     THE MEMORY REQUIREMENTS JUST BELOW DON'T MATCH THE ACTUAL CODE.
C
C---      MXNEED= NEEDA+NEEDB+NANGM*NOC2*NBF
C---      MINEED= NEEDA+NEEDB+NOC2*NBF
C---      DMXNEED=(NOC2*NVIR)/125000.0D+00
C---      IF(GOPARR) THEN
C---C        MXDISK= REAL(NOC2)*REAL(NVIR)*REAL(MXP)/125000.0
C---         MXDISK= INT(DMXNEED*MXP)
C---      ELSE
C---C        MXDISK= REAL(NOC2)*REAL(NVIR)*REAL(NBF)/125000.0
C---         MXDISK= INT(DMXNEED*NBF)
C---      ENDIF
C---      IF(MASWRK) WRITE(IW,8) MXNEED,MINEED,MXDISK
C---    8 FORMAT(/1X,'MEMORY AND DISK REQUIREMENTS FOR SEMI-
C---     *        'DIRECT MP2 TRANSFORMATION'/
C---     *        '  MAXIMUM MEMORY=',I10,' WORDS'/
C---     *        '  MINIMUM MEMORY=',I10,' WORDS'/
C---     *        '  REQUIRED DISK =',I10,' MBYTES')
C
C        TOTAL DISK SPACE IS NOC * NOC * NVIR * NBF, AVOID OVERFLOWING
C
      DMXNEED = NOC2*NVIR
      DMXNEED = DMXNEED/125000.0D+00
      MXDISK  = INT(NBF*DMXNEED)
      MXDISKP = INT(NBF*DMXNEED/NPROC)
      IF(MASWRK  .AND.  ICHKPASS.EQ.1)
     *     WRITE(IW,9020) MXDISK,MXDISKP,NPROC
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
C     MAXC = IDAMAX(NBF*NQMT,X(IVEC),1)
      MAXC = IDAMAX(NBF*NO,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     ----- ASSIGN -IRREP- LABEL TO MO'S -----
C
      CALL VALFMP(LOADFM)
      IWRK  = 1     + LOADFM
      IWRK2 = IWRK  + NBF3
      IWRK3 = IWRK2 + NBF2
      LAST  = IWRK3 + NBF
      NEEDC = LAST -LOADFM - 1
      CALL GETFMP(NEEDC)
      IF(ICHKPASS.EQ.1) GO TO 100
 
      CALL DAREAD(IDAF,IODA,X(IWRK),NBF3,12,0)
      CALL DAREAD(IDAF,IODA,X(IWRK2),NBF2,45,0)
      CALL TRFSYM(X(ILAB),X(IIRP),X(IDEG),X(IWRK2),X(IWRK),
     *            X(IVEC),X(IWRK3),IA,NO,NBF,NO,NBF)
C
  100 CONTINUE
      CALL RETFMP(NEEDC)
C
C     ----- TRANSPOSE MO COEFFICIENTS AND REORDER VIRTUAL MOS ------
C
      CALL MOTRN(X(IVEC),X(IVECO),X(IVECV),X(IIRP),X(INSYM),X(IISYM),
     *           NOC,NVIR)
C
C     ----- ALLOCATE MEMORY FOR 2ND TRANSFORMATION -----
C
      CALL GOTFMP(NGOTMX)
      JPTR2 = NANGM*NOC2*NBF
      DIVANG=.FALSE.
      IF(NGOTMX.LT.JPTR2) THEN
         DO 110 I=1, NANGM-1
            MXANGM= NANGM-I
            JPTR2= MXANGM*NOC2*NBF
            IF(NGOTMX.GE.JPTR2) GOTO 120
  110    CONTINUE
         IF(MASWRK) WRITE(IW,9040) NEEDA+JPTR2
         CALL ABRT
  120    DIVANG=.TRUE.
         IF(MASWRK) WRITE(IW,9050) MXANGM
      ENDIF
      CALL VALFMP(LOADFM)
      IPTR2   = 1       + LOADFM
      LBSMAX  = IPTR2   + JPTR2
      LLBLINT = LBSMAX  + NBF
      LNUMINT = LLBLINT + NBF*NANGM**3
      LAST    = LNUMINT + NANGM**3
      NEEDD = LAST - LOADFM -1 
C
      NEED123 = NEEDA + NEEDB + NEEDD
      IF(MASWRK) THEN
         IF(ICHKPASS.EQ.1) WRITE(IW,9060) NEED123
         IF(ICHKPASS.EQ.2) WRITE(IW,9065)
         CALL FLSHBF(6)
      END IF
C
      CALL GETFMP(NEEDD)
      IF(EXETYP.EQ.CHECK) GO TO 130
C
C     ----- OPEN DIRECT ACCESS FILE(S) ON EACH NODE -----
C
      IDAF20= 20
      IF(MASWRK) CALL GMS_GETENV('DASORT',PATHNM)
      IF(ISGDDI) THEN
        PATHNM=ENVBUF(IDAF20)
        IF(PAROUT.AND.NODEXT(IDAF20).EQ.0) CALL ADDNANODE(PATHNM,MEGLOB)
      ELSE
         CALL PARENV('DASORT',PATHNM,IOUT)
      ENDIF
      NULL = CHAR(0)
      DO 3 KOL=1,256
         IF(PATHNM(KOL:KOL).EQ.' '  .OR.
     *      PATHNM(KOL:KOL).EQ.NULL) GO TO 4
    3 CONTINUE
      KOL=257
    4 CONTINUE
      IF(KOL.EQ.1) THEN
         WRITE(IW,1) 'DASORT'
         CALL ABRT
    1 FORMAT(1X,'YOU MUST ASSIGN GENERIC NAME ',A,' WITH A SETENV.')
      END IF
      KOL=KOL-1
      OPEN(UNIT=IDAF20,FILE=PATHNM(1:KOL),STATUS='UNKNOWN',
     *     ACCESS='DIRECT', FORM='UNFORMATTED', RECL=8*NVIR)
C
C     ----- CARRY OUT 1ST, 2ND, AND 3RD INDEX TRANSFORMATIONS -----
C
      CALL MP2SD1(X(IVECO),X(IVECV),X(IGHND),
     *            X(IXINT),X(IINT),X(IPTR1),X(IPTR2),X(IPTR3),X(IWORK),
     *            X(LBSMAX),X(LLBLINT),X(LNUMINT),
     *            NUMP,NOC,NVIR,NOC2,NANGM,CUTOFF,
     *            NLOOP,ILOOP,DIVANG,MXANGM,NIJ3)
C
  130 CONTINUE
      CALL RETFMP(NEEDD)
      CALL RETFMP(NEEDB)
C
C     ----- ALLOCATE MEMORY FOR 4TH TRANSFORMATION -----
C
      CALL GOTFMP(NGOTMX)
      NGOTMX = NGOTMX - 2*NVIR*NVIR - NOC2 - NVIR*NBF
      IF(GOPARR) THEN
         IF(ICHKPASS.EQ.1) THEN
            MXP=NBF/NPROC
         ELSE
         IF(DLB) THEN
            MXP=0
            DO I= 1,NPROC
               MXP= MAX(MXP,NLOOP(2,I))
            ENDDO
         ENDIF
         END IF
         MIJ4= NOC2/NPROC
         IF(MOD(NOC2,NPROC).NE.0) MIJ4= MIJ4+1
         NIJ4 = NGOTMX/(NVIR*(NBF+MXP*NPROC))
         IF(NIJ4.GT.MIJ4) NIJ4= MIJ4
         IF(NIJ4.EQ.0) THEN
           IF(MASWRK) WRITE(IW,*) 'PARALLEL MP2IMS: NOT ENOUGH MEMORY'
           CALL ABRT
         ENDIF
         JDDIS= NIJ4*NVIR*MXP*NPROC
         NPASS= NOC2/(NIJ4*NPROC)
         IF(MOD(NOC2,(NIJ4*NPROC)).NE.0) NPASS= NPASS+1
      ELSE
         NIJ4 = NGOTMX/(NVIR*NBF)
         IF(NIJ4.GT.NOC2) NIJ4= NOC2
         IF(NIJ4.EQ.0) THEN
           WRITE(IW,*) 'SERIAL MP2IMS: NOT ENOUGH MEMORY'
           CALL ABRT
         ENDIF
         JDDIS= 1
         NPASS= NOC2/NIJ4
         IF(MOD(NOC2,NIJ4).NE.0) NPASS= NPASS+1     
      ENDIF
C
      CALL VALFMP(LOADFM)
      IDDIS = 1     + LOADFM
      JPTR3 = IDDIS + JDDIS
      IMO1  = JPTR3 + NIJ4*NVIR*NBF
      IENGO = IMO1  + NVIR*NVIR
      IWORK = IENGO + NOC2
      LAST  = IWORK + NVIR*NBF
      NEEDE = LAST  - LOADFM -1
C
      NEED4 = NEEDA + NEEDE
      IF(MASWRK  .AND.  ICHKPASS.EQ.1) WRITE(IW,9070) NEED4
C
      CALL GETFMP(NEEDE)
      IF(EXETYP.EQ.CHECK) GO TO 230
C
C     ----- START 4TH QUARTER TRANSFORMATION -----
C     -----       AND MP2 ENERGY CALCULATION -----
C
C
      CALL MOCHNG(X(IVECV),X(IVECV2),NVIR,NLOOP,ILOOP)
      CALL MP2QT4(X(IENG),X(IENGO),X(IVECV2),X(IWORK),X(IDDIS),
     *            X(JPTR3),X(IMO1),X(IIRP),X(INSYM),X(IISYM),
     *            NVIR,NOC,NOC2,NT,
     *            IJEND,ILEN,IRECV,NLOOP,NCPU,MXP,NPASS,NIJ4,E2S,E2T)
C
C     ----- DELETE THE DIRECT ACCESS FILE(S) ON EACH CORE -----
C
      CLOSE(UNIT=IDAF20, STATUS='DELETE')
C
  230 CONTINUE
      CALL RETFMP(NEEDE)
      CALL RETFMP(NEEDA)
C
C     ----- PRINT ENERGY DATA -----
C
      IF(GOPARR) CALL DDI_GSUMF(1000,E2S,1)
      IF(GOPARR) CALL DDI_GSUMF(1000,E2T,1)
      IF(ICHKPASS.EQ.2  .OR.  EXESAV.EQ.CHECK) THEN
         EHF =ESCF
         E2=E2S+E2T
         EMP2=EHF+E2
         SCSE2= E2S*P12+E2T/THREE
         ESCSMP2=EHF+SCSE2
         IF(MASWRK) WRITE(IW,9080) EHF,E2,EMP2,E2S,E2T,SCSE2,ESCSMP2
         IF(SCSPT.NE.ANONE) THEN
            E2=SCSE2 
            EMP2=ESCSMP2
         ENDIF
      END IF
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
      IF(.NOT.ABEL) CALL SYMON
      RETURN
C
 9000 FORMAT(/1X,26('-'),3X,46(1H-)/
     *        1X,'RHF-MP2 ENERGY CALCULATION',3X,
     *           'PROGRAM WRITTEN BY K.ISHIMURA'/
     *        1X,26('-'),3X,46(1H-)/
     *        1X,'NUMBER OF OCCUPIED ORBITALS            =',I5,/
     *        1X,'NUMBER OF CORRELATED OCCUPIED ORBITALS =',I5,/
     *        1X,'NUMBER OF VIRTUAL ORBITALS             =',I5,/
     *        1X,'NUMBER OF BASIS FUNCTIONS              =',I5,/
     *        1X,'NUMBER OF BASIS SHELLS                 =',I5,/
     *        1X,'INTEGRAL THRESHOLD                     =',1P,E10.2)
 9010 FORMAT(/1X,'NUMBER OF REDUNDANT MOS IN FMO         =',I5)
 9020 FORMAT(1X,'TOTAL DISK REQUIRED (ALL PROCESSORS)=',I15,' GBYTES'/
     *       1X,'                 DISK SPACE PER CORE=',I15,
     *          ' GBYTES, USING P=',I5)
 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 RHF-MP2 *****',/,
     *       1X,'INCREASE MEMORY TO AT LEAST',I10,'WORDS.')
 9050 FORMAT(1X,'MULTIPLE PASSES',
     *      /1X,'MAXIMUM ANGULAR MOMENTUM=',I3)
 9060 FORMAT(1X,'MEMORY REQUIRED FOR TRANSFORMATION OF 1ST THREE',
     *          ' INDICES=',I12)
 9065 FORMAT(1X,'BEGINNING MP2 INTEGRAL TRANSFORMATION...')
 9070 FORMAT(1X,'MEMORY REQUIRED FOR  FINAL INDEX TRANSF. AND MP2',
     *          ' ENERGY=',I12)
 9080 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//,
     *       13X,' E(2S)=',1X,F20.10/,
     *       13X,' E(2T)=',1X,F20.10/,
     *       13X,'E(2ST)=',1X,F20.10/,
     *       13X,'SCSMP2=',1X,F20.10)
      END
C
C*MODULE MP2IMS  *DECK MOTRN
      SUBROUTINE MOTRN(VEC,TRNOCC,TRNVIR,MOSIRP,NSYM,ISYM,NOC,NVIR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
      DIMENSION VEC(NBF,*),TRNOCC(NOC,*),TRNVIR(NBF,*),MOSIRP(*),
     *          NSYM(*),ISYM(NVIR,*)
C
C     TRANSPOSE OCCUPIED AND VIRTURAL MO COEFFICIENTS 
C
      DO 110 I=1,NBF
         DO 110 J=1,NOC
            TRNOCC(J,I)= VEC(I,J+NACORE)
  110 CONTINUE
C
      IF(NT.NE.1) THEN
         DO 130 IT=1,NT
            ICOUNT=0
            DO 120 MO=1,NVIR
               IF(MOSIRP(MO+NOA).EQ.IT) THEN
                 ICOUNT=ICOUNT+1
                 ISYM(ICOUNT,IT)=MO
               ENDIF
  120       CONTINUE
            NSYM(IT)=ICOUNT
  130    CONTINUE
         NMO=1
         DO 150 IT=1,NT
            DO 150 J=1,NSYM(IT)
               MO= ISYM(J,IT)
               DO 140 I=1,NBF
                  TRNVIR(I,NMO)= VEC(I,MO+NOA)
  140          CONTINUE
               NMO=NMO+1
  150    CONTINUE
      ELSE
        DO 160 I=1,NVIR
           DO 160 J=1,NBF
              TRNVIR(J,I)= VEC(J,I+NOA)
  160   CONTINUE
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2IMS  *DECK MOCHNG
      SUBROUTINE MOCHNG(VECV,VECV2,NVIR,NLOOP,ILOOP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      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 VECV(NBF,*),VECV2(NBF,*)
      DIMENSION NLOOP(2,NPROC),ILOOP(NSHELL,NPROC)
C
C     EXCHANGE VIRTUAL MO COEFFICIENTS BEFORE 4TH TRANSFORMATION
C
      IF(GOPARR.AND.IBTYP.EQ.1) THEN
         N=0
         DO 110 I=1,NPROC
            DO 110 J=1,NLOOP(1,I)
               ISH= ILOOP(J,I)
               MINI= KMIN(ISH)
               MAXI= KMAX(ISH)
               LOCI= KLOC(ISH)-MINI
               DO 110 K=MINI,MAXI
                  N=N+1
                  DO 110 L=1,NVIR
                     VECV2(N,L)=VECV(K+LOCI,L)
  110    CONTINUE
      ELSE
         CALL DCOPY(NVIR*NBF,VECV,1,VECV2,1)
      ENDIF
      RETURN
      END
C
C
C*MODULE MP2IMS  *DECK MP2SD1
      SUBROUTINE MP2SD1(VECO,VECV,GHONDO,
     *                  XINTS,AOINTS,PTR1,PTR2,PTR3,WORK,
     *                  BSMAX,LBLINT,NUMINT,
     *                  NUMP,NOC,NVIR,NOC2,NANGM,CUTOFF,
     *                  NLOOP,ILOOP,DIVANG,MXANGM,NIJ3)
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 /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
      COMMON /SHLNOS/ QQ4,IDUMMY(20)
C
      DIMENSION VECO(NOC,*),VECV(*),
     *          GHONDO(*),XINTS(*),AOINTS(*),PTR1(*),PTR2(*),
     *          PTR3(*),WORK(*),NUMP(0:NPROC)
      DIMENSION NLOOP(2,0:NPROC-1),ILOOP(NSHELL,0:NPROC-1)
      DIMENSION BSMAX(NBF),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
      DO I=1,NBF
         DUM=ZER
         DO J=1,NOC
            DUM=MAX(DUM,ABS(VECO(J,I)))
         ENDDO
         BSMAX(I)=DUM
      ENDDO
      JSTART=1
      DO I=1,NSHELL
         DUM=ZER
         DO J= JSTART, KLOC(I)
            DUM=MAX(DUM,BSMAX(J))
         ENDDO
         JSTART= KLOC(I)+1
      ENDDO
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(1,ME)= 0
         NLOOP(2,ME)= 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(1,ME)= NLOOP(1,ME)+1
            ILOOP(NLOOP(1,ME),ME)=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(2,ME)=NLOOP(2,ME)+NUMI
         CALL VCLR(PTR2,1,NOC2*NBF*NUMI)
C
C     ----- START K SHELL -----
C
      DO 400 KSH= 1,NSHELL
         CALL SYMIK(ISH,KSH,Q2,NSHELL)
         IF(Q2.EQ.ZER) GO TO 400
         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
         CALL SYMIKLJ(ISH,KSH,LSH,JSH,Q4,NSHELL)
         IF(Q4.EQ.ZER) GO TO 200
         QQ4=Q4
         IJIJ= IA(ISH)+JSH
         IF(JSH.GT.ISH) IJIJ= IA(JSH)+ISH
         TEST= QQ4*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 MP2QT1(VECO,AOINTS,PTR1,NOC,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 MP2QT2(VECO,PTR1,PTR2,
     *            NOC,NOC2,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 MP2QT3(VECV,PTR2,PTR3,WORK,
     *            NVIR,NOC2,IREC,NUMI,NIJ3)
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
C     ----- CHECK DYNAMIC PARALLEL INFORMATIONS -----
      IF(NXT.AND.GOPARR) THEN
         CALL DDI_DLBRESET
         INTSIZ=8/NWDVAR
C        INTSIZ=4
CI32     INTSIZ=4
CI64     INTSIZ=8
         DO 600 I=1,NPROC-1
C           IF(MASWRK .AND. I.NE.MASTER) THEN
            IF(MASWRK) THEN
               CALL DDI_RECV(NLOOP(1,I),2*INTSIZ,I)
               CALL DDI_RECV(ILOOP(1,I),NSHELL*INTSIZ,I)
            ELSE
               IF(I.EQ.ME)CALL DDI_SEND(NLOOP(1,I),2*INTSIZ,MASTER)
               IF(I.EQ.ME)CALL DDI_SEND(ILOOP(1,I),NSHELL*INTSIZ,MASTER)
            ENDIF
  600    CONTINUE
         CALL DDI_BCAST(1012,'I',NLOOP,2*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,'TH SHELL LOOP START')
      RETURN
      END
C
C
C*MODULE MP2IMS  *DECK SAVEAO
      SUBROUTINE SAVEAO(AOINTS,GHONDO,CUTOFF,
     *                  LBLINT,NUMINT,KSH,LSH,LOCJ,MINI,
     *                  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 /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      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
      IKL  = 0
      LMAX = MAXL
C
      DO 240 I= MINI2,MAXI2
         I_INDEX = (I-MINI)*LSTRI + 1
         DO 230 K= MINK,MAXK
            IK_INDEX = (K-MINK)*LSTRK + I_INDEX
            IF(KANDL) LMAX= K
            DO 220 L= MINL,LMAX
               IKL_INDEX = (L-MINL)*LSTRL + IK_INDEX
               IKL= IKL+1
               DO 210 J= MINJ,MAXJ
                  JJ= LOCJ+J
                  IJKL_INDEX = (J-MINJ)*LSTRJ + IKL_INDEX
C
                  VAL = GHONDO( IJKL_INDEX)
CISHI             IF(ABS(VAL)*BSMAX(JJ).LT.CUTOFF) GO TO 210
                  IF(ABS(VAL).LT.CUTOFF) GO TO 210
                  IF(KANDL.AND.K.EQ.L) VAL= VAL*HALF
                  NUMINT(IKL)= NUMINT(IKL)+1
                  AOINTS(NUMINT(IKL),IKL)= VAL
                  LBLINT(NUMINT(IKL),IKL)= JJ
  210          CONTINUE
  220       CONTINUE
  230    CONTINUE
  240 CONTINUE
C
      RETURN
      END
C
C
C*MODULE MP2IMS  *DECK MP2QT1
      SUBROUTINE MP2QT1(VECO,AOINTS,PTR1,NOC,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
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
C
      DIMENSION VECO(NOC,*),AOINTS(NBF,*),
     *          PTR1(NOC,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,NOC
                  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,NOC
                     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*MODULE MP2IMS  *DECK MP2QT2
      SUBROUTINE MP2QT2(VECO,PTR1,PTR2,
     *                  NOC,NOC2,NUMI,NUMK,NUML,NANGM,CUTOFF,
     *                  KSH,LSH,MINK,MINL,LOCK,LOCL)

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL KNOTL
C
      PARAMETER (SQRT2=1.414213562373095D+00)
C
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
C
      DIMENSION VECO(NOC,*),PTR1(NOC,NANGM,NANGM,NANGM),PTR2(NOC2,NBF,*)
C
C     2ND QUARTER TRANSFORMATION
C
      KNOTL= KSH.NE.LSH
      LOCK2= LOCK+MINK-1
      LOCL2= LOCL+MINL-1
      
C
      IF(KNOTL)THEN
      DO 160 I= 1,NUMI
         DO 140 L= 1,NUML
            LL= LOCL2+L
            DO 150 K= 1,NUMK
               KK= LOCK2+K
               MOIJ= 0
               DO 130 MOI= 1,NOC
                  DUM1= PTR1(MOI,L,K,I)
CISHI             IF(ABS(DUM1)*TMP1.LT.CUTOFF) GO TO 120
                  IF(ABS(DUM1).LT.CUTOFF) GO TO 120
                  DUM2= DUM1*SQRT2
                  DO 110 MOJ= 1,MOI-1
                     PTR2(MOIJ+MOJ,LL,I)= PTR2(MOIJ+MOJ,LL,I)
     *                                   +DUM2*VECO(MOJ,KK)
  110             CONTINUE
                  PTR2(MOIJ+MOI,LL,I)= PTR2(MOIJ+MOI,LL,I)
     *                                +DUM1*VECO(MOI,KK)
  120             MOIJ= MOIJ+MOI
  130          CONTINUE
  150       CONTINUE
  140    CONTINUE
  160 CONTINUE
C
      DO 260 I=1,NUMI
         DO 250 K=1,NUMK
            KK= LOCK2+K
            DO 240 L=1,NUML
               LL= LOCL2+L
               MOIJ= 0
               DO 230 MOI=1,NOC
                  DUM1= PTR1(MOI,L,K,I)
CISHI             IF(ABS(DUM1)*TMP1.LT.CUTOFF) GO TO 220
                  IF(ABS(DUM1).LT.CUTOFF) GO TO 220
                  DUM2= DUM1*SQRT2
                  DO 210 MOJ= 1,MOI-1
                     PTR2(MOIJ+MOJ,KK,I)= PTR2(MOIJ+MOJ,KK,I)
     *                                   +DUM2*VECO(MOJ,LL)
  210             CONTINUE
                  PTR2(MOIJ+MOI,KK,I)= PTR2(MOIJ+MOI,KK,I)
     *                                +DUM1*VECO(MOI,LL)
  220             MOIJ= MOIJ+MOI
  230          CONTINUE
  240       CONTINUE
  250    CONTINUE
  260 CONTINUE
C
      ELSE
      DO 380 I=1,NUMI
         DO 370 K= 1,NUMK
            LMAX= K
            KK= LOCK2+K
            DO 360 L= 1,LMAX
               LL= LOCL2+L
               MOIJ= 0
               DO 350 MOI= 1,NOC
                  DUM1= PTR1(MOI,L,K,I)
                  DUM2= DUM1*SQRT2
                  DUM3= ABS(DUM1)
CISHI             IF(DUM3*TMP1.LT.CUTOFF) GO TO 320
                  IF(DUM3.LT.CUTOFF) GO TO 320
                  DO 310 MOJ= 1,MOI-1
                     PTR2(MOIJ+MOJ,LL,I)= PTR2(MOIJ+MOJ,LL,I)
     *                                   +DUM2*VECO(MOJ,KK)
  310             CONTINUE
                  PTR2(MOIJ+MOI,LL,I)= PTR2(MOIJ+MOI,LL,I)
     *                                +DUM1*VECO(MOI,KK)
CI320             IF(DUM3*TMP2.LT.CUTOFF) GO TO 340
  320             IF(DUM3.LT.CUTOFF) GO TO 340
                  DO 330 MOJ= 1,MOI-1
                     PTR2(MOIJ+MOJ,KK,I)= PTR2(MOIJ+MOJ,KK,I)
     *                                   +DUM2*VECO(MOJ,LL)
  330             CONTINUE
                  PTR2(MOIJ+MOI,KK,I)= PTR2(MOIJ+MOI,KK,I)
     *                                +DUM1*VECO(MOI,LL)
  340             MOIJ= MOIJ+MOI
  350          CONTINUE
  360       CONTINUE
  370    CONTINUE
  380 CONTINUE
      ENDIF
C
      RETURN
      END
C
C
C*MODULE MP2IMS  *DECK MP2QT3
      SUBROUTINE MP2QT3(VECV,PTR2,PTR3,WORK,
     *                  NVIR,NOC2,IREC,NUMI,NIJ3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER(ZER=0.0D+00, ONE=1.0D+00)
C
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /MP2DAF/ IDAF20,NAV20
C
      DIMENSION VECV(NBF,*),PTR2(NOC2,NBF,*),PTR3(NVIR,NIJ3),
     *          WORK(NBF,*)
C
C     3RD QUARTER TRANSFORMATION
C
      IPASS= NOC2/NIJ3
      MODIJ= MOD(NOC2,NIJ3)
      IF(MODIJ.NE.0) THEN
         IPASS= IPASS+1
      ELSE
         MODIJ= NIJ3
      ENDIF
C
      DO 150 I= 1,NUMI
         MIJ3= NIJ3
         DO 140 IJPASS= 1,IPASS
            IF(IJPASS.EQ.IPASS) MIJ3= MODIJ
            MOIJ2= (IJPASS-1)*NIJ3
            DO 120 K= 1,NBF
               DO 110 MOIJ= 1,MIJ3
                  WORK(K,MOIJ)= PTR2(MOIJ+MOIJ2,K,I)
  110          CONTINUE
  120       CONTINUE
            CALL DGEMM('T','N',NVIR,MIJ3,NBF,ONE,VECV,NBF,
     *                 WORK,NBF,ZER,PTR3,NVIR)
            DO 130 IJWRIT= 1,MIJ3
               IREC= IREC+1
               WRITE(UNIT=IDAF20,REC=IREC)
     *              (PTR3(MOB,IJWRIT),MOB=1,NVIR)
  130       CONTINUE
  140    CONTINUE
  150 CONTINUE
C
      RETURN
      END
C
C
C*MODULE MP2IMS  *DECK MP2QT4
      SUBROUTINE MP2QT4(ENG,ENGO,VECV,WORK,DDIS,
     *                  PTR3,PTR4,MOSIRP,NSYM,ISYM,NVIR,NOC,NOC2,NT,
     *                  IJEND,ILEN,IRECV,NLOOP,NCPU,
     *                  MXP,NPASS,NIJ4,E2S,E2T)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SYM2EI
C
      PARAMETER (MXAO=8192)
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
C
C        NOTE THAT -NCPU- IS AN ALIAS FOR -NPROC-
      DIMENSION IJEND(NCPU),ILEN(NCPU),IRECV(NCPU),NLOOP(2,0:NCPU-1)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /MP2DAF/ IDAF20,NAV20
C
      DIMENSION ENG(NBF),ENGO(NOC2),VECV(NBF,*),
     *          DDIS(NVIR,NIJ4,MXP,*),PTR3(NVIR,NIJ4,*),PTR4(NVIR,*),
     *          WORK(NBF,*),MOSIRP(*),NSYM(*),ISYM(NVIR,*)
C
C     4TH QUARTER TRANSFORMATION AND MP2 ENERGY CALCULATION
C
C     ----- CALCULATE I+J MO ENERGIES -----
C
      DO 120 I= 1,NOC
         ENGI= ENG(I+NACORE)
         II  = IA(I)
         DO 110 J= 1,I
            ENGO(II+J)= ENGI+ENG(J+NACORE)
  110    CONTINUE
  120 CONTINUE
C
C     ----- 4TH QUARTER TRANSFORMATION -----
C     ----- AND MP2 ENERGY CALCULATION -----
C
C     ----- PARALLEL VERSION -----
C
      IF(GOPARR) THEN
         NUMI=NLOOP(2,ME)
         LENDDI= NVIR*NIJ4*NUMI
         MIJ4  = NIJ4
         MODIJ1= MOD(NOC2,NIJ4*NPROC)
         MODIJ2= MOD(MODIJ1,NPROC)
         MIJTMP= MODIJ1/NPROC
         IF(MIJTMP.EQ.0.AND.MODIJ2.EQ.0) MIJTMP= NIJ4
C
         NSUM=1
         DO 210 IPROC= 1,NPROC
            IJEND(IPROC)= MIJTMP
            IF(IPROC.LE.MODIJ2) IJEND(IPROC)= IJEND(IPROC)+1
            ILEN(IPROC)=NLOOP(2,IPROC-1)*NVIR*NIJ4
            IRECV(IPROC)=NSUM
            NSUM= NSUM+NLOOP(2,IPROC-1)
  210    CONTINUE
C
C     WRITE(6,*) ME,'MP2 CODE REACHED THE DATA EXCHANGE.' 
      DO 500 IPASS= 1,NPASS
         IJREC1=(IPASS-1)*NIJ4*NPROC
         MIJ4= NIJ4
         DO 330 I= 1,NUMI
            IJREC2= (I-1)*NOC2+IJREC1
            DO 320 IPROC=1,NPROC
               IF(IPASS.EQ.NPASS) MIJ4= IJEND(IPROC)
               DO 310 MOIJ= 1,MIJ4
                 READ(UNIT=IDAF20,REC=IJREC2+MOIJ)
     *                (DDIS(MOB,MOIJ,I,IPROC),MOB=1,NVIR)
  310          CONTINUE
               IJREC2= IJREC2+MIJ4
  320       CONTINUE
  330    CONTINUE
C
C     ----- CALL NON-BLOCKING SEND AND RECIEVE ---
C
         JPROC= ME+1
         DO 340 IPROC= ME+2,NPROC
            JPROC= JPROC-1
            IF(JPROC.EQ.0) JPROC= NPROC
            CALL DDI_ISEND(DDIS(1,1,1,IPROC),LENDDI*8,IPROC-1,
     *                        IREQ1)
            CALL DDI_IRECV(PTR3(1,1,IRECV(JPROC)),ILEN(JPROC)*8,
     *                     JPROC-1,IREQ2)
            CALL DDI_WAIT(IREQ2)
            CALL DDI_WAIT(IREQ1)
  340    CONTINUE
         DO 350 IPROC= 1,ME
            JPROC= JPROC-1
            IF(JPROC.EQ.0) JPROC= NPROC
            CALL DDI_ISEND(DDIS(1,1,1,IPROC),LENDDI*8,IPROC-1,
     *                        IREQ1)
            CALL DDI_IRECV(PTR3(1,1,IRECV(JPROC)),ILEN(JPROC)*8,
     *                     JPROC-1,IREQ2)
            CALL DDI_WAIT(IREQ2)
            CALL DDI_WAIT(IREQ1)
  350    CONTINUE
         CALL DCOPY(LENDDI,DDIS(1,1,1,ME+1),1,PTR3(1,1,IRECV(ME+1)),1)
C
         MIJ4= NIJ4
         MIJ5= ME*NIJ4
         IF(IPASS.EQ.NPASS) THEN
            MIJ4= IJEND(ME+1)
            MIJ5= 0
            DO 360 IPROC= 1,ME
               MIJ5= MIJ5+ IJEND(IPROC)
  360       CONTINUE
         ENDIF
         IJREC2=MIJ5+(IPASS-1)*NIJ4*NPROC
C
         IF(NT.NE.1) THEN
            DO 370 II= 1,NOC-1
               MOI= II
               IF(IJREC2.LT.IA(II+1)) GO TO 380
  370       CONTINUE
            MOI= NOC
  380       MOJ= IJREC2-IA(MOI)+1
         ENDIF
C
         DO 490 MOIJ=1,MIJ4
            DO 400 I= 1,NBF
               DO 390 MOB= 1,NVIR
                  WORK(I,MOB)= PTR3(MOB,MOIJ,I)
  390          CONTINUE
  400       CONTINUE
            EIJ= ENGO(IJREC2+MOIJ)
C
            IF(NT.EQ.1) THEN
               CALL DGEMM('T','N',NVIR,NVIR,NBF,ONE,WORK,NBF,
     *                    VECV,NBF,ZER,PTR4,NVIR)
               DO 420 MOA= 1,NVIR
                  EMO1= ENG(MOA+NOA)-EIJ
                  DO 410 MOB= 1,NVIR
                     XAB= PTR4(MOB,MOA)
                     IF(ABS(XAB).LT.TOL) GO TO 410
                     XBA= PTR4(MOA,MOB)
                     EMO2= ONE/(EMO1+ENG(MOB+NOA))
                     E2S= E2S-XAB*XAB*EMO2
                     E2T= E2T-XAB*(XAB-XBA)*EMO2
  410             CONTINUE
  420          CONTINUE
C
            ELSE
               NIT= 1
               DO 440 IT= 1,NT
                  IF(NSYM(IT).EQ.0) GO TO 440
                  NJT= 1
                  DO 430 JT= 1,NT
                     IF(NSYM(JT).EQ.0) GO TO 430
                     MOSA= ISYM(1,IT)+NOA
                     MOSB= ISYM(1,JT)+NOA
                     IF(SYM2EI(MOSIRP(MOI+NACORE),MOSIRP(MOSA),
     *                         MOSIRP(MOJ+NACORE),MOSIRP(MOSB)))
     *                  CALL DGEMM('T','N',NSYM(JT),NSYM(IT),NBF,ONE,
     *                             WORK(1,NJT),NBF,VECV(1,NIT),NBF,ZER,
     *                             PTR4(1,NIT),NVIR)
  430             NJT= NJT+NSYM(JT)
  440          NIT= NIT+NSYM(IT)
               NIT= 0
               DO 480 IT= 1,NT
                  IF(NSYM(IT).EQ.0) GO TO 480
                  NJT= 0
                  DO 470 JT= 1,NT
                     IF(NSYM(JT).EQ.0) GO TO 470
                     MOSA= ISYM(1,IT)+NOA
                     MOSB= ISYM(1,JT)+NOA
                     IF(SYM2EI(MOSIRP(MOI+NACORE),MOSIRP(MOSA),
     *                         MOSIRP(MOJ+NACORE),MOSIRP(MOSB))) THEN
                        DO 460 MOA= 1,NSYM(IT)
                           EMO1= ENG(ISYM(MOA,IT)+NOA)-EIJ
                           DO 450 MOB= 1,NSYM(JT)
                              XAB= PTR4(MOB,MOA+NIT)
                              IF(ABS(XAB).LT.TOL) GO TO 450
                              XBA= PTR4(MOA,MOB+NJT)
                              EMO2= ONE/(EMO1+ENG(ISYM(MOB,JT)+NOA))
                              E2S= E2S-XAB*XAB*EMO2
                              E2T= E2T-XAB*(XAB-XBA)*EMO2
  450                      CONTINUE
  460                   CONTINUE
                     ENDIF
  470             NJT= NJT+NSYM(JT)
  480          NIT= NIT+NSYM(IT)
               MOJ= MOJ+1
               IF(MOI.LT.MOJ) THEN
                  MOI= MOI+1
                  MOJ= 1
               ENDIF
            ENDIF
  490    CONTINUE
  500 CONTINUE
C
C     ----- SERIAL VERSION -----
C
      ELSE
      MIJ4 = NIJ4
      DO 660 IPASS= 1,NPASS
         IJPASS=(IPASS-1)*NIJ4
         IF(IPASS.EQ.NPASS.AND.(MOD(NOC2,NIJ4).NE.0))
     *      MIJ4= MOD(NOC2,NIJ4)
         DO 520 I= 1,NBF
            IJREC= (I-1)*NOC2+(IPASS-1)*NIJ4
            DO 510 MOIJ= 1,MIJ4
               READ(UNIT=IDAF20,REC=IJREC+MOIJ)
     *             (PTR3(MOB,MOIJ,I),MOB=1,NVIR)
  510       CONTINUE
  520    CONTINUE

         IF(NT.NE.1) THEN
            DO 530 II= 1,NOC-1
               MOI=II
               IF(IJPASS.LT.IA(II+1)) GO TO 540
  530       CONTINUE
            MOI= NOC
  540       MOJ= IJPASS-IA(MOI)+1
         ENDIF
C
         DO 650 MOIJ= 1,MIJ4
            DO 560 I= 1,NBF
               DO 550 MOB= 1,NVIR
                  WORK(I,MOB)= PTR3(MOB,MOIJ,I)
  550          CONTINUE
  560       CONTINUE
            EIJ= ENGO(MOIJ+IJPASS)
C
            IF(NT.EQ.1) THEN
               CALL DGEMM('T','N',NVIR,NVIR,NBF,ONE,WORK,NBF,
     *                    VECV,NBF,ZER,PTR4,NVIR)
               DO 580 MOA= 1,NVIR
                  EMO1= ENG(MOA+NOA)-EIJ
                  DO 570 MOB= 1,NVIR
                     XAB= PTR4(MOB,MOA)
                     IF(ABS(XAB).LT.TOL) GO TO 570
                     XBA = PTR4(MOA,MOB)
                     EMO2= ONE/(EMO1+ENG(MOB+NOA))
                     E2S = E2S-XAB*XAB*EMO2
                     E2T = E2T-XAB*(XAB-XBA)*EMO2
  570             CONTINUE
  580          CONTINUE
C
            ELSE
               NIT=1
               DO 600 IT=1,NT
                  IF(NSYM(IT).EQ.0) GO TO 600
                  NJT=1
                  DO 590 JT=1,NT
                     IF(NSYM(JT).EQ.0) GO TO 590
                     MOSA= ISYM(1,IT)+NOA
                     MOSB= ISYM(1,JT)+NOA
                     IF(SYM2EI(MOSIRP(MOI+NACORE),MOSIRP(MOSA),
     *                         MOSIRP(MOJ+NACORE),MOSIRP(MOSB)))
     *                  CALL DGEMM('T','N',NSYM(JT),NSYM(IT),NBF,ONE,
     *                             WORK(1,NJT),NBF,VECV(1,NIT),NBF,ZER,
     *                             PTR4(1,NIT),NVIR)
  590             NJT= NJT+NSYM(JT)
  600          NIT= NIT+NSYM(IT)
               NIT= 0
               DO 640 IT= 1,NT
                  IF(NSYM(IT).EQ.0) GO TO 640
                  NJT= 0
                  DO 630 JT= 1,NT
                     IF(NSYM(JT).EQ.0) GO TO 630
                     MOSA= ISYM(1,IT)+NOA
                     MOSB= ISYM(1,JT)+NOA
                     IF(SYM2EI(MOSIRP(MOI+NACORE),MOSIRP(MOSA),
     *                         MOSIRP(MOJ+NACORE),MOSIRP(MOSB))) THEN
                        DO 620 MOA= 1,NSYM(IT)
                           EMO1= ENG(ISYM(MOA,IT)+NOA)-EIJ
                           DO 610 MOB= 1,NSYM(JT)
C                             XAB= PTR4(MOB,MOA)
                              XAB= PTR4(MOB,MOA+NIT)
                              IF(ABS(XAB).LT.TOL) GO TO 610
C                             XBA= PTR4(MOA,MOB)
                              XBA= PTR4(MOA,MOB+NJT)
                              EMO2= ONE/(EMO1+ENG(ISYM(MOB,JT)+NOA))
                              E2S= E2S-XAB*XAB*EMO2
                              E2T= E2T-XAB*(XAB-XBA)*EMO2
  610                      CONTINUE
  620                   CONTINUE
                     ENDIF
  630             NJT= NJT+NSYM(JT)
  640          NIT= NIT+NSYM(IT)
               MOJ= MOJ+1
               IF(MOI.LT.MOJ) THEN
                  MOI= MOI+1
                  MOJ= 1
               ENDIF
            ENDIF
  650    CONTINUE
  660 CONTINUE
      ENDIF
C
      RETURN
      END
C
C*MODULE MP2IMS  *DECK SHELLQUARTMP2
      SUBROUTINE SHELLQUARTMP2(GHONDO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GHONDO(*)
      DIMENSION GROTSPD(1296)
C
      LOGICAL PACK2E,IANDJ,KANDL,SAME
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH, MXSH=5000, MXGTOT=20000)
C
      COMMON /FLIPS / IB(4,3)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      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 /SHLEXC/ NORGSH(3),NORGSP(3),IEXCH,NANGM,NGTH(4)
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
C
C  FOR RYS QUADRATURE CODE
C
      LOGICAL LSHEL
      COMMON /FMCOM /XX(1)
C
C  FOR ERIC CODE, BUT ALL CODES MUST SET ERIOUT FOR OUTPUT ROUTINES
C
      LOGICAL ERICQT
      COMMON /ERIDAT/ LEN1,LEN2,LEN3,LEN4
      COMMON /ERIOUT/ INW,JNW,KNW,LNW,LSTRI,LSTRJ,LSTRK,LSTRL
C
C  FOR THE TWO ROTATED AXIS CODES
C
      INTEGER IDPOP(4,10)
      LOGICAL SP,SPD,ROTSP,ROTSPD
      COMMON /GOUT  / GPOPLE(768),NORGP
      COMMON /POPOUT/ LPOPI,LPOPJ,LPOPK,LPOPL
      DATA IDPOP/0,0,0,0,216,36,6,1,432,72,12,2,648,108,18,3,
     *           0,0,0,0,216,36,6,1,432,72,12,2,648,108,18,3,
     *           864,144,24,4,1080,180,30,5/
C
C     ----- SELECT THE INTEGRAL CODE FOR THIS SHELL QUARTET -----
C     THE USER INPUT SELECTION -INTTYP- HAS THE FOLLOWING MEANING:
C      INTTYP=0    BEST TIMING:
C                  USE ROTATED AXIS CODES FOR ANY S,P,D,L SHELLS,
C                  OTHERWISE PICK ERIC CODE WHENEVER POSSIBLE, BUT
C                  USE RYS CODE FOR L SHELL OR IF TOO MUCH ANG.MOM.
C      INTTYP=1    USE THE S,P,L OR S,P,D,L ROTATED AXIS CODE WHENEVER
C                  POSSIBLE, OTHERWISE RYS QUADRATURE (NO ERIC).
C      INTTYP=2    USE ERIC CODE AS MUCH AS POSSIBLE, OTHERWISE
C                  USE THE RYS QUADRATURE (NO ROTATED AXIS).
C      INTTYP=3    USE RYS POLYNOMIAL QUADRATURE FOR EVERYTHING.
C
      SP    = KTYPE(ISH).LE.2.AND.
     *        KTYPE(JSH).LE.2.AND.
     *        KTYPE(KSH).LE.2.AND.
     *        KTYPE(LSH).LE.2
      SPD   = KTYPE(ISH).LE.3.AND.
     *        KTYPE(JSH).LE.3.AND.
     *        KTYPE(KSH).LE.3.AND.
     *        KTYPE(LSH).LE.3
      IF(SP) SPD=.FALSE.
      LSHEL = (KMAX(ISH)-KMIN(ISH)+1).EQ.4.OR.
     *        (KMAX(JSH)-KMIN(JSH)+1).EQ.4.OR.
     *        (KMAX(KSH)-KMIN(KSH)+1).EQ.4.OR.
     *        (KMAX(LSH)-KMIN(LSH)+1).EQ.4
      LQSUM = KTYPE(ISH) + KTYPE(JSH) + KTYPE(KSH) + KTYPE(LSH) - 4
C
      ROTSP  = SP
      ROTSPD = SPD
      ERICQT = .NOT.LSHEL  .AND.  LQSUM.LE.5
C
C        RYS QUADRATURE IS A BIT FASTER AT UNCONTRACTED QUARTETS
C
      KQCON = KNG(ISH) * KNG(JSH) * KNG(KSH) * KNG(LSH)
      IF(KQCON.EQ.1) THEN
         ROTSP  = .FALSE.
         ROTSPD = .FALSE.
         ERICQT = .FALSE.
      END IF
C
C        INPUT OVERRIDES
C
      IF(INTTYP.EQ.1) ERICQT = .FALSE.
      IF(INTTYP.GE.2) ROTSP  = .FALSE.
      IF(INTTYP.GE.2) ROTSPD = .FALSE.
      IF(INTTYP.EQ.3) ERICQT = .FALSE.
C     THE CODE BELOW KNOWS NOT HOW TO HANDLE ERIC.
      ERICQT = .FALSE.
C
C        THE VARIOUS PACKAGES ARE CALLED BELOW IN THE ORDER OF
C        FIRST ROTATED AXIS, THEN ERIC, FINALLY RYS QUADRATURE.
C        SINCE EACH CODE RETURNS, A QUARTET IS NEVER DONE TWICE.
C
C  ROTATED AXIS CODE FOR PURE SP SHELL QUARTET
C
      IF (ROTSP) THEN
        INW = ISH
        JNW = JSH
        KNW = KSH
        LNW = LSH
        NORGP = NORGSP(IEXCH)
        NORGH = NORGSH(IEXCH)
C
        CALL GENR70(1,.FALSE.)
C
C  SAVE TO OUTPUT ARRAY WITH HONDO INDEXING
C
        MINI = KMIN(INW)
        MAXI = KMAX(INW)
        MINJ = KMIN(JNW)
        MAXJ = KMAX(JNW)
        MINK = KMIN(KNW)
        MAXK = KMAX(KNW)
        MINL = KMIN(LNW)
        MAXL = KMAX(LNW)
C
        II = 1
        DO I = MINI, MAXI
          IP = (I-1)*LPOPI + 1
          IJ  = II
          DO J = MINJ, MAXJ
            IJP = (J-1)*LPOPJ + IP
            IJK  = IJ
            DO K = MINK, MAXK
              IJKP = (K-1)*LPOPK + IJP
              IJKL  = IJK
              DO L = MINL, MAXL
                IJKLP = (L-1)*LPOPL + IJKP
                GHONDO(IJKL+NORGH) = GPOPLE(IJKLP+NORGP)
                IJKL = IJKL  + LEN1
              END DO
              IJK  = IJK  + LEN2
            END DO
            IJ  = IJ  + LEN3
          END DO
          II = II + LEN4
        END DO
        LSTRI = LEN4
        LSTRJ = LEN3
        LSTRK = LEN2
        LSTRL = LEN1
        RETURN
C       ******
C
C  ROTATED AXIS CODE FOR QUARTET CONTAINING AT LEAST ONE D FUNCTION
C
      ELSE IF (ROTSPD) THEN
        INW = ISH
        JNW = JSH
        KNW = KSH
        LNW = LSH
C
        CALL GENR03(GROTSPD)
C
C  SAVE TO OUTPUT ARRAY WITH HONDO INDEXING
C
        NORGH = NORGSH(IEXCH)
        IANDJ = ISH.EQ.JSH
        KANDL = KSH.EQ.LSH
        SAME  = ISH.EQ.KSH  .AND.  JSH.EQ.LSH
        IF(NOPK.EQ.0) SAME=.FALSE.
C
        IEX=1
        IBB = IB(1,IEX)
        JBB = IB(2,IEX)
        KBB = IB(3,IEX)
        LBB = IB(4,IEX)
C
        MINI = KMIN(INW)
        MAXI = KMAX(INW)
        MINJ = KMIN(JNW)
        MAXJ = KMAX(JNW)
        MINK = KMIN(KNW)
        MAXK = KMAX(KNW)
        MINL = KMIN(LNW)
        MAXL = KMAX(LNW)
C
        IJN = 0
        JMAX = MAXJ
        DO I = MINI, MAXI
          IHONDO = (I-MINI)*LEN4 + 1
          IROTAX = IDPOP(IBB,I)  + 1
          DO 340 J = MINJ, JMAX
            IJHONDO = (J-MINJ)*LEN3 + IHONDO
            IJROTAX = IDPOP(JBB,J)  + IROTAX
            IJN = IJN+1
            LMAX=MAXL
            KLN=0
            DO K = MINK, MAXK
              IJKHONDO = (K-MINK)*LEN2 + IJHONDO
              IJKROTAX = IDPOP(KBB,K)  + IJROTAX
              DO L = MINL, LMAX
                KLN = KLN+1
                IJKLHONDO = (L-MINL)*LEN1 + IJKHONDO
                IJKLROTAX = IDPOP(LBB,L)  + IJKROTAX
                GHONDO(IJKLHONDO+NORGH) = GROTSPD(IJKLROTAX)
              END DO
            END DO
  340     CONTINUE
        END DO
        LSTRI = LEN4
        LSTRJ = LEN3
        LSTRK = LEN2
        LSTRL = LEN1
        RETURN
C       ******
C
C  USE ERIC FAST CODES, REQUIRES THAT LQSYM.LE.5 AND NO L-SHELLS
C
      ELSE IF (ERICQT) THEN
        NORGH = NORGSH(IEXCH)
        CALL ERIC(ISH,JSH,KSH,LSH,GHONDO(1+NORGH))
        RETURN
C       ******
C
C  GENERAL CASE = HONDO/RYS QUADRATURE: ANY S,P,D,F,G, OR L SHELLS
C
      ELSE
        CALL VALFM(LOADFM)
        IDDIJ = LOADFM + 1
C       NEED  = 16*MXG2
        NEED  = 49*MXG2
        CALL GETFM(NEED)
        CALL SHELLS(1,ISH,JSH,KSH,LSH,.FALSE.)
        CALL SHELLS(2,ISH,JSH,KSH,LSH,.FALSE.)
        CALL IJPRIM(XX(IDDIJ))
        NORGH = NORGSH(IEXCH)
        IF(NOPK.NE.0) CALL ZQOUT(GHONDO)
        IF(LQSUM.EQ.0) THEN
           CALL S0000(GHONDO(1+NORGH),XX(IDDIJ))
        ELSE
           CALL GENRAL(GHONDO(1+NORGH),XX(IDDIJ))
        END IF
        CALL RETFM(NEED)
C
        IANDJ = ISH.EQ.JSH
        KANDL = KSH.EQ.LSH
        SAME  = ISH.EQ.KSH  .AND.  JSH.EQ.LSH
        IF(NOPK.EQ.0) SAME=.FALSE.
C
        IF(IANDJ.OR.SAME) THEN
          MINI = KMIN(ISH)
          MAXI = KMAX(ISH)
          MINJ = KMIN(JSH)
          MAXJ = KMAX(JSH)
          MINK = KMIN(KSH)
          MAXK = KMAX(KSH)
          MINL = KMIN(LSH)
          MAXL = KMAX(LSH)
          IJN = 0
          JMAX = MAXJ
          DO I = MINI, MAXI
            IHND = I-MINI
            IF(IANDJ) JMAX=I
            DO J = MINJ, JMAX
              JHND = J-MINJ 
              IJN = IJN+1
              LMAX=MAXL
              KLN=0
              DO K = MINK, MAXK
                KHND = K-MINK
                IF(KANDL) LMAX=K
                DO L = MINL, LMAX
                  KLN = KLN+1
                  LHND = L-MINL
                  IJKLHND = IHND*LEN4+JHND*LEN3+KHND*LEN2+LHND*LEN1+1
                  IF(SAME .AND. IJN.GT.KLN) THEN
                    IF(IANDJ) THEN
                      JIKLHND= JHND*LEN4+IHND*LEN3+KHND*LEN2+LHND*LEN1+1
                      KLIJHND= KHND*LEN4+LHND*LEN3+IHND*LEN2+JHND*LEN1+1
                      LKIJHND= LHND*LEN4+KHND*LEN3+IHND*LEN2+JHND*LEN1+1
                      GHONDO(JIKLHND+NORGH) = GHONDO(IJKLHND+NORGH)
                      GHONDO(KLIJHND+NORGH) = GHONDO(IJKLHND+NORGH) 
                      GHONDO(LKIJHND+NORGH) = GHONDO(IJKLHND+NORGH) 
                    ELSE
                      KLIJHND= KHND*LEN4+LHND*LEN3+IHND*LEN2+JHND*LEN1+1
                      GHONDO(KLIJHND+NORGH) = GHONDO(IJKLHND+NORGH)
                    ENDIF
                  ELSEIF(IANDJ) THEN
                    JIKLHND = JHND*LEN4+IHND*LEN3+KHND*LEN2+LHND*LEN1+1
                    GHONDO(JIKLHND+NORGH) = GHONDO(IJKLHND+NORGH) 
                  ENDIF
                END DO
              END DO
            END DO
          END DO
        ENDIF
        RETURN
C       ******
C
      END IF
      END
C
C*MODULE UNPORT  *DECK SETFMP
      SUBROUTINE SETFMP(IBUF,NBUF)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /FMPARP/ LTOP,LOFFS,MEMLIM
C
C     ----- SET THE PRIVATE FM ARRAY -----
C     IBUF SHOULD POINT TO A BUFFER ALREADY ALLOCATED BY GETFM.
C     (I.E., IBUF IS RETURNED BY LOADFM, WITH THE FURTHER
C     ADDITION OF 1 AND, POSSIBLY, OTHER OFF-SETS). 
C     IN OTHER WORDS, SETFMP PERFORMES MEMORY ALLOCATION WITHIN AN 
C     ALREADY ALLOCATED MEMORY. IF NONE WAS OBTAINED (IBUF=0), THEN
C     THE STANDARD FUNCTIONS WIILL BE CALLED. 
C
      LTOP=0
      LOFFS=IBUF
      MEMLIM=NBUF
      RETURN
      END
C
C*MODULE UNPORT  *DECK VALFMP
      SUBROUTINE VALFMP(IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /FMPARP/ LTOP,LOFFS,MEMLIM
C
C     ----- RETURN THE CURRENT TOP OF FM ARRAY -----
C
      IF(LOFFS.NE.0) THEN 
         IPAR = LTOP + LOFFS
C        WRITE(6,*) 'VALFMP: ',LTOP,LOFFS,IPAR 
      ELSE  
         CALL VALFM(IPAR)
      ENDIF  
      RETURN
      END
C*MODULE UNPORT  *DECK GETFMP
      SUBROUTINE GETFMP(IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /FMPARP/ LTOP,LOFFS,MEMLIM
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- RESERVE IPAR WORDS OF FM -----
C
      IF(LOFFS.NE.0) THEN
         LTOP = LTOP + IPAR
         IF (LTOP .GT. MEMLIM ) THEN
            WRITE(IW,9000) ME,LTOP,MEMLIM
            CALL ABRT
C           ONE COULD ALLOCATE FROM GETFM; BUT THEN ONE ALSO
C           SHOULD REMEMBER THAT WHEN RETURNING...
         ENDIF
C        WRITE(6,*) 'GETFMP: ',LTOP,IPAR
      ELSE
         CALL GETFM(IPAR)
      ENDIF
      RETURN
C
 9000 FORMAT(1X,'***** ERROR: PRIVATE MEMORY REQUEST EXCEEDS AVAILABLE',
     *          ' MEMORY'/,
     *     1X,'PROCESS NO.',I5,' WORDS REQUIRED=',I10,' AVAILABLE=',I10)
      END
C*MODULE UNPORT  *DECK RETFMP
      SUBROUTINE RETFMP(IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /FMPARP/ LTOP,LOFFS,MEMLIM
C
C     ----- RETURN IPAR WORDS OF FM -----
C
      IF(LOFFS.NE.0) THEN
         LTOP = LTOP - IPAR
         IF(LTOP.LT.0) THEN
            WRITE(6,*) 'MEMORY PROBLEM, RETFMP: ',LTOP,IPAR
            CALL ABRT
         END IF
      ELSE
         CALL RETFM(IPAR)
      ENDIF
      RETURN
      END
C*MODULE UNPORT  *DECK GOTFMP
      SUBROUTINE GOTFMP(IPAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /FMPARP/ LTOP,LOFFS,MEMLIM
C
C     ----- RETURN NUMBER OF FREE WORDS IN FM -----
C
      IF(LOFFS.NE.0) THEN
         IPAR = MEMLIM - LTOP
C        WRITE(6,*) 'GOTFMP: ',IPAR,MEMLIM,LTOP
      ELSE
         CALL GOTFM(IPAR)
      ENDIF
      RETURN
      END
C
C*MODULE UNPORT  *DECK SETFMPI
      SUBROUTINE SETFMPI(MINMEM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK,PACK2E
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /PCKLAB/ LABSIZ
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     USE THE INTEGRAL BUFFER AS THE PRIVATE FAST MEMORY.
C     THERE IS NO NEED TO "UNSET" IT. CAVEAT EMPTOR TEMENS!
C     DO NOT USE THE BUFFER IF IT IS SMALL ( < MEMPRI IN $MP2).
C
      L2EBUF=0
      N2EBUF=0 
      IF(NINTIC.GT.0) THEN
        LABSIZ2=2/LABSIZ
        N2EBUF=NINTMX+NINTIC+(NINTMX+NINTIC-1)/LABSIZ2+1
        IF(N2EBUF.GE.MINMEM.AND.MINMEM.GT.0) THEN 
          L2EBUF=LBUFPIC
          IF(MASWRK) WRITE(IW,9000) N2EBUF
        ELSE
          N2EBUF=0
        ENDIF
      ENDIF
      CALL SETFMP(L2EBUF,N2EBUF)
C     WRITE(6,*) 'SETFMPI: ',LBUFPIC,MINMEM,L2EBUF,N2EBUF
      RETURN
 9000 FORMAT(1X,'USING ',I9,' WORDS OF THE INCORE 2E INTEGRAL BUFFER.')
      END
C
C*MODULE MP2NEW  *DECK SYMIK
      SUBROUTINE SYMIK(II,KK,Q2,NSHELL)
C
C -----------------------------------------------------------------
C  DETERMINE SYMMETRY-UNIQUE SHELL PAIRS
C -----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000)
      PARAMETER (ZER=0.0D+00, P12=1.0D-12)
C
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
      Q2= ZER
      IK=(II-1)*NSHELL+KK
C
C  LOOP THRU OPERATIONS IN THE GROUP AND MAP TO NEW PAIRS
C
      N= 1
      DO IT= 2,NT
C
C  MAP CENTERS
C
         II_NEW= MAPSHL(II,IT)
         KK_NEW= MAPSHL(KK,IT)
C
C  COMPARE CANONICAL PAIR INDICES
C
         IK_NEW=(II_NEW-1)*NSHELL+KK_NEW
         IF(IK.LT.IK_NEW) RETURN
         IF(IK.EQ.IK_NEW) N= N+1
      ENDDO
      Q2= NT
      Q2= Q2/N
      IF(ABS(Q2-NINT(Q2)).GT.P12) CALL ABRT
C
      RETURN
      END
C
C*MODULE MP2NEW  *DECK SYMIKLJ
      SUBROUTINE SYMIKLJ(II,KK,LL,JJ,Q4,NSHELL)
C
C -----------------------------------------------------------------
C  DETERMINE SYMMETRY-UNIQUE QUARTETS
C -----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000,MXATM=2000, MXAO=8192)
      PARAMETER (ZER=0.0D+00, P12=1.0D-12)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
      Q4= ZER
      NS2= NSHELL*(NSHELL*NSHELL+NSHELL)/2
C
C  ORIGINAL 4-INDEX IN THE NON-CANONICAL LIST KK.GE.LL
C
      INDX= (II-1)*NS2+(IA(KK)+LL-1)*NSHELL+JJ
C
C  LOOP THRU OPERATIONS IN THE GROUP AND MAP TO NEW PAIRS
C
      N= 1
      DO IT= 2,NT
C
C  MAP CENTERS
C
         II_NEW = MAPSHL(II,IT)
         JJ_NEW = MAPSHL(JJ,IT)
         KK_NEW = MAPSHL(KK,IT)
         LL_NEW = MAPSHL(LL,IT)
C
C  MAPPED 4-INDEX IN THE NON-CANONICAL LIST KK.GE.LL
C
         IF(KK_NEW.GE.LL_NEW) THEN
            INDX_NEW=(II_NEW-1)*NS2+(IA(KK_NEW)+LL_NEW-1)*NSHELL+JJ_NEW
         ELSE
            INDX_NEW=(II_NEW-1)*NS2+(IA(LL_NEW)+KK_NEW-1)*NSHELL+JJ_NEW
         ENDIF
C
C  COMPARE INDICES
C
         IF(INDX_NEW.GT.INDX) RETURN
         IF(INDX_NEW.EQ.INDX) N= N+1
      ENDDO
      Q4 = NT
      Q4 = Q4/N
      IF(ABS(Q4-NINT(Q4)).GT.P12) CALL ABRT
C
      RETURN
      END
