C 22 DEC 06 - DGF - SYNCHRONISE MP2PAR
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 13 MAR 06 - MWS - SEPARATE DAF RECORD NOS OF RHF/ROHF MP2 FROM UHF MP2
C 22 FEB 06 - MWS - IMPLEMENT OPTION FOR NO SYMMETRY IN MP2 STEP
C  2 FEB 06 - MWS - PROPERTY CALLS MOVED TO WFNMP2 DRIVER
C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - MWS - ADJUSTMENT FOR H OR I FUNCTION RUNS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - SYNCHRONISE SYMBLK COMMON BLOCK
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - SPW - AOCPHF,FLMAT,FLMATU: CHANGES FOR PARALLEL CIS
C  3 SEP 03 - CMA - CHANGES TO ALLOW CPHF ITERATION LIMIT INPUT
C 12 DEC 02 - CMA - MP2GE2,MP2GC2,MP2GRD,UMPGE2,UMPGC2,UMPGRD:CPHF TWEAK
C  7 AUG 02 - MWS - MP2GRD: FIX MEMORY LEAK FOR COSMO
C 22 MAY 02 - CMA - LOOP REORDERINGS FOR EFFICIENCY
C 29 APR 02 - MWS - AOCPHF,CPHFU: ALLOW TO CONTINUE IF CONVERGING
C 24 JAN 02 - CMA - ADD SOME ROUTINES FOR UMP2 GRADIENTS
C 13 JUN 01 - MWS - SPARSQ: REMOVE PAUSE STATEMENTS, MPCPHF: FIX MROT=0
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 11 JUN 00 - MWS - MP2GE2,FLMAT: ALWAYS ALLOCATE HONDO INTEGRAL BUFFER
C 10 APR 00 - MWS - REMOVE STATIC MEMORY FROM COSMO COMMONS
C 25 MAR 00 - KKB/LNB - MP2 COSMO PROPERTIES/OCE CORRECTION
C 16 FEB 00 - MWS - MP2GRD: FIX MPPROP MEMORY LEAK
C 21 DEC 99 - MWS - CHANGE INTEGRAL SYMMETRY TEST
C  9 APR 99 - MWS - MP2GRD: PRINT MESSAGES BEFORE ABRT, UPDATE /MP2PAR/
C 13 FEB 99 - MWS - AOCPHF: SUPPRESS ITERATION PRINTOUT
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C  3 NOV 98 - MWS - MP2GRD: COPY MP2 ENERGY TO ENRGYS COMMON
C 27 OCT 98 - MAF - MAMNLS,MP2GE2,MP2GRD: ALLOW USE OF SPHER. HARMONICS
C  6 JAN 98 - MWS - CHANGE RACLOS CALL
C 24 JUN 97 - MWS - AOCPHF: REMOVE UNUSED ARGUMENT FROM CPDIIS CALL
C 27 MAR 97 - MWS - MP2GRD: RAISE ITERMX=30 UP TO 50
C 18 DEC 96 - TLW - AOCPHF, CPDIIS: MODIFICATIONS TO AO DRIVEN CPHF
C  9 DEC 96 - JHJ - AOCPHF: NEW GUESS, FORCE 2 DIIS ITERATIONS
C 29 SEP 96 - MWS - FLMAT,MP2GE2: USE ONLY NECESSARY GHONDO BUFFER SIZE
C 11 JUL 96 - MWS - MP2GRD: AVOID 2PDM FOR PROPERTY RUNS, GENERATE NO-S
C 27 JUN 96 - TLW - AOCPHF: RETURN FAST MEMORY
C 13 JUN 96 - SPW - ADDED AOCPHF,CPDIIS,FLMAT,W2IJAO, PLUS OTHER CHANGES
C                   FOR FROZEN CORE MP2 GRADIENTS IN MP2GE2,MP2GC2,
C                   P2W2N,MP2GRD,MP2LAI,MAMNLS
C  9 FEB 96 - MWS - MP2GRD: FIX MPCPHO MEMORY OVERALLOCATION BUG
C 11 NOV 95 - NM,SPW,MWS - FINISH IMPLEMENTATION OF HONDO MP2 GRADIENTS
C
C*MODULE MP2GRD  *DECK JKDMP2
      SUBROUTINE JKDMP2(PHF,PMP2,C,XIN,NBF3,NBF2,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION PHF(NBF2),PMP2(NBF2),XIN(NBF2,NBF,*),C(NBF3)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2DM2/ NOC1,NOC2,ISTEP,NOC3,NOC4
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA NFT16/16/
      DATA RHF,ROHF/8HRHF     ,8HROHF    /
C
C        READ PARTIALLY TRANSFORMED -DM2- FROM DISK
C
      IF(SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.ROHF) THEN
         KREC7 = 307
         KREC8 = 308
      ELSE
         KREC7 = 417
         KREC8 = 418
      END IF
C
      IF(ISTEP.EQ.1) CALL SEQREW(NFT16)
      DO 15 MI=1,NOC2-NOC1+1
         DO 10 NMU=1,NBF
            CALL TRFRD(NFT16,XIN(1,NMU,MI),NX,PHF,NIX,NBF2)
   10    CONTINUE
   15 CONTINUE
C
C        READ SCF AND MP2 FIRST ORDER DENSITY MATRICES
C
      IF(ISTEP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,PHF ,NBF2,KREC8,0)
         CALL DAREAD(IDAF,IODA,PMP2,NBF2,KREC7,0)
         CALL DAREAD(IDAF,IODA,C   ,NBF3,   15,0)
      END IF
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MAMNLS
      SUBROUTINE MAMNLS(AIBJ,E,C,X,Y,XX,XIN,NBF,NBF2,NOC,
     *                  NVIR,NFTA,NFT16,TOL,SOME,NCO,NOCA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
      PARAMETER (TWO=2.0D+00)
C
      LOGICAL SOME
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /IJPAIR/ IA(MXAO)
C
      DIMENSION AIBJ(NVIR,NVIR,NOCA),C(NBF,*),X(NBF,*),Y(NBF,*),
     *          XIN(NBF2,NBF),XX(*),E(*)
C
C         GENERATE NON-SEPARABLE CONTRIBUTIONS TO MP2 2ND ORDER DENSITY
C
      IF(SOME) WRITE(IW,9000)
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTA)
      NO = 0
C
      DO 1000 MI=1,NOCA
C
C     ----- READ ALL (AI|BJ) INTEGRALS FOR GIVEN -I- -----
C
      DO 115 MA=1,NVIR
         CALL TRFRD(NFTA,XX,NX,X,NIX,NBF2)
         DO 110 MB=1,NVIR
            MBJ=IA(MB+NOC)+NCO
            DO 105 MJ=1,NOCA
               AIBJ(MB,MA,MJ) = XX(MBJ+MJ)
  105       CONTINUE
  110    CONTINUE
  115 CONTINUE
C
      CALL VCLR(XIN,1,NBF*NBF2)
C
C     ----- -X(A,B)- IS -GAMMA(A,I,B,J)- FOR GIVEN -I,J- -----
C
      DO 150 MJ=1, NOCA
         DO 135 MB=1,NVIR
            DO 130 MA=1,NVIR
               DIJAB = E(MI+NCO) + E(MJ+NCO)
     *                           - E(MA+NOC) - E(MB+NOC)
               X(MA,MB) = TWO*(TWO*AIBJ(MB,MA,MJ)-
     *                    AIBJ(MA,MB,MJ))/DIJAB
  130       CONTINUE
  135    CONTINUE
C
C     ----- TRANSFORM -A,B- TO -MU,LAMBDA- -----
C
         CALL MRARTR(X,NBF,NVIR,NVIR,C(1,NOC+1),NBF,NBF,Y,NBF)
         CALL MRARBR(C(1,NOC+1),NBF,NBF,NVIR,Y,NBF,NBF,X,NBF)
C
C     ----- TRANSFORM -J- TO -SIGMA- AND SYMMETRIZE -SIGMA,LAMBDA -----
C     ----- -XIN(MU,LASI)- IS -GAMMA(MU,I,LASI)- FOR GIVEN -I-    -----
C
         DO 148 NMU=1, NBF
            LASI = 0
            DO 143 NLA=1, NBF
               DO 140 NSI=1, NLA
                  LASI = LASI+1
                  XIN(LASI,NMU)=XIN(LASI,NMU) + C(NSI,MJ+NCO)*X(NMU,NLA)
     *                                        + C(NLA,MJ+NCO)*X(NMU,NSI)
  140          CONTINUE
  143       CONTINUE
  148   CONTINUE
  150 CONTINUE
C
      DO 160 NMU=1, NBF
         CALL TRFWTM(NFT16,XIN(1,NMU),NX,X,NIX,NBF2,1,TOL)
         NO = NO+NX
  160 CONTINUE
C
 1000 CONTINUE
C
      IF(SOME) WRITE(IW,9999) NO,NFT16
      RETURN
C
 9000 FORMAT(/1X,'FORMING NON-SEPARABLE PART OF THE MP2 2ND ORDER',
     *           ' DENSITY MATRIX...')
 9999 FORMAT(1X,I12,' -2DM- ELEMENTS WRITTEN ON FILE',I3)
      END
C
C*MODULE MP2GRD  *DECK MPCPHF
      SUBROUTINE MPCPHF(WXY,WXY2,XX,IX,E,MOSIRP,XLAI,NORB,NORB2,
     *                  NOC,NVIR,T,B,NROT,NFTI,NFTA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,THREE=3.0D+00,FOUR=4.0D+00)
C
      LOGICAL SOME,OUT,DBUG
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
C
      EXTERNAL SPARAX1,SPARXA1
C
      DIMENSION WXY(NOC,NVIR,NOC,NVIR),WXY2(NVIR*NOC,NVIR*NOC)
      DIMENSION E(NORB),MOSIRP(NORB),XLAI(NVIR,NOC)
      DIMENSION XX(NORB2),IX(NORB2)
      DIMENSION T(NROT,4),B(NROT)
C
C     ----- IN-CORE SOLUTION TO THE CPHF EQUATIONS OF MP2 -----
C
      CALL SEQREW(NFTI)
      CALL SEQREW(NFTA)
C
      DO 15 I=1,NROT
         DO 10 J=1,NROT
            WXY2(J,I) = ZERO
   10    CONTINUE
   15 CONTINUE
C
C ... (IJ||AB) INTEGRALS
C
      DO 45 J=1,NOC
         DO 40 I=J,NOC
            CALL TRFRD(NFTI,XX,NX,IX,NIX,NORB2)
            IF(DBUG) CALL TRFPRT(XX,IX,NORB,I,J)
            DO 35 K=1,NVIR
               KL = IA(K+NOC)+NOC
               DO 20 L=1,K
                  KL=KL+1
                  WXY(I,K,J,L) = WXY(I,K,J,L) - XX(KL)
   20          CONTINUE
               IF(I.NE.J) THEN
                  KL = IA(K+NOC)+NOC
                  DO 30 L=1,K-1
                     KL=KL+1
                     WXY(J,K,I,L) = WXY(J,K,I,L) - XX(KL)
   30             CONTINUE
               END IF
   35       CONTINUE
   40    CONTINUE
   45 CONTINUE
C
C ... (AI||BJ) INTEGRALS
C
      DO 65 J=1,NOC
         DO 60 I=1,NVIR
            CALL TRFRD(NFTA,XX,NX,IX,NIX,NORB2)
            IF(DBUG) CALL TRFPRT(XX,IX,NORB,I+NOC,J)
            DO 55 K=1,I-1
               KL = IA(K+NOC)
               DO 50 L=1,NOC
                  KL=KL+1
                  WXY(J,I,L,K) = WXY(J,I,L,K) + FOUR*XX(KL)
                  WXY(L,I,J,K) = WXY(L,I,J,K) - XX(KL)
   50          CONTINUE
   55       CONTINUE
            KL = IA(I+NOC)
            DO 58 L=1,J
               KL=KL+1
               WXY(J,I,L,I) = WXY(J,I,L,I) + THREE*XX(KL)
   58       CONTINUE
   60    CONTINUE
   65 CONTINUE
C
C     ----- ENERGY TERM TO DIAGONAL -----
C
      DO 85 I=1,NVIR
         DO 80 J=1,NOC
            WXY(J,I,J,I) = WXY(J,I,J,I) + E(I+NOC) - E(J)
   80    CONTINUE
   85 CONTINUE
C
C     ----- COPY LOWER-LEFT PART IN UPPER-RIGHT -----
C
      DO 95 IJ=1,NROT-1
         DO 90 KL=IJ+1,NROT
            WXY2(IJ,KL) = WXY2(KL,IJ)
   90    CONTINUE
   95 CONTINUE
C
      IF(DBUG) THEN
         WRITE(IW,9999)
         CALL PRSQ(WXY2,NROT,NROT,NROT)
      END IF
C
C ... KEEP ONLY NON-ZERO ELEMENTS
C
      IF(NIRRED.EQ.1) THEN
         MROT = NROT
         IJ=0
         DO 105 I=1,NVIR
            DO 100 J=1,NOC
               IJ=IJ+1
               B(IJ) = XLAI(I,J)
  100       CONTINUE
  105    CONTINUE
      ELSE
         MROT = 0
         IJ   = 0
         DO 135 I=1,NVIR
            DO 130 J=1,NOC
               IJ=IJ+1
               IF(MOSIRP(I+NOC).EQ.MOSIRP(J)) THEN
                  MROT = MROT+1
                  DO 110 KL=1,NROT
                     WXY2(MROT,KL) = WXY2(IJ,KL)
  110             CONTINUE
                  DO 120 KL=1,NROT
                     WXY2(KL,MROT) = WXY2(KL,IJ)
  120             CONTINUE
                  B(MROT) = XLAI(I,J)
               END IF
  130       CONTINUE
  135    CONTINUE
C
         IF(DBUG) THEN
            WRITE(IW,9998)
            CALL PRSQ(WXY,MROT,MROT,NROT)
         END IF
      END IF
C
C       AN EXTREMELY UNLIKELY BRANCH, E.G. H2 IN MBS, WHERE THE
C       VIRTUAL IS OF DIFFERENT SYMMETRY THAN THE OCCUPIED ORBITAL
C
      IF(MROT.EQ.0) GO TO 165
C
      DO 140 I=1,MROT
         XX(I) = ZERO
         T(I,4) = ONE/WXY2(I,I)
  140 CONTINUE
C
      DO 160 I=1,MROT
         DO 150 J=1,MROT
            WXY2(J,I) = WXY2(J,I)*T(J,4)
  150    CONTINUE
         B(I) = B(I)*T(I,4)
  160 CONTINUE
C
      IF(SOME) WRITE(IW,9996) MROT
      CALL SPARSQ(WXY,NROT,B,XX,MROT,T(1,1),T(1,2),T(1,3),T(1,4),
     *            SPARAX1,SPARXA1)
C
  165 CONTINUE
      IF(NIRRED.EQ.1) THEN
         IJ=0
         DO 175 I=1,NVIR
            DO 170 J=1,NOC
               IJ=IJ+1
               XLAI(I,J) = XX(IJ)
  170       CONTINUE
  175    CONTINUE
      ELSE
         IJ=0
         KL=0
         DO 185 I=1,NVIR
            DO 180 J=1,NOC
               IJ=IJ+1
               IF(MOSIRP(I+NOC).EQ.MOSIRP(J)) THEN
                  KL=KL+1
                  XLAI(I,J) = XX(KL)
               ELSE
                  XLAI(I,J) = ZERO
               END IF
  180       CONTINUE
  185    CONTINUE
      END IF
C
      IF(DBUG) THEN
         WRITE(IW,9997)
         CALL PRSQ(XLAI,NOC,NVIR,NVIR)
      END IF
C
      RETURN
C
 9999 FORMAT(/' MP2-CPHF MATRIX (FULL)')
 9998 FORMAT(/' MP2-CPHF MATRIX (REDUCED)')
 9997 FORMAT(/' PAI(VIRT,OCC)')
 9996 FORMAT(5X,'CHOOSING IN CORE SOLUTION OF THE CPHF EQUATIONS,'/
     *       5X,'THE SIZE OF THE LINEAR EQUATION SYSTEM IS',I10)
      END
C
C*MODULE MP2GRD  *DECK MPCPHO
      SUBROUTINE MPCPHO(WXY,IIJJ,XX,IX,E,MOSIRP,XLAI,NORB,
     *                  NORB2,NOC,NVIR,T,B,NROT,NFTI,NFTA,CUTOFF,MEM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,THREE=3.0D+00,FOUR=4.0D+00)
C
      LOGICAL SOME,OUT,DBUG
C
C     NOTE THAT -IIJJ- AND -B- ARE EQUIVALENT BY CALL.
C     -IX- AND -XX- ARE USED WITH LENGTHS OF -NORB2- AND -NROT- HERE.
C
      DIMENSION WXY(MEM),IIJJ(NROT),IX(*),XX(*),E(NORB),
     *          XLAI(NVIR,NOC),B(NROT),T(NROT,5),MOSIRP(NORB)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /MP2PTR/ NFT,J00,J01,J02,IDUMMY(4)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
C
      EXTERNAL MPSPAX,MPOUAX
C
C     ----- OUT-OF-CORE SOLUTION TO THE CPHF EQUATIONS FOR MP2 -----
C     ROWS OF THE ORBITAL HESSIAN WILL BE WRITTEN TO FILE -NFT-
C
      NFT  = 16
      CALL SEQREW(NFT)
C
C     ----- LOOP UNTIL ALL -WXY- ELEMENTS HAVE BEEN CALCULATED -----
C
      NO   = 0
      MROT = 0
      MINIJ= 1
 1000 CONTINUE
C
C     ----- SET UP THIS SLICE OF THE ORBITAL HESSIAN -----
C     THE HESSIAN -WXY- IS COMPUTED USING TRIANGULAR STORAGE.
C     COMPLETE ROWS FROM -MINIJ- TO -MAXIJ- ARE INCLUDED.
C     THE SET OF ROWS COMPUTED MUST FIT WITHIN -MEM- WORDS.
C
      IJ=0
      KL=MINIJ
   10 CONTINUE
      IIJJ(KL)=IJ
      IJ=IJ+KL
      IF(IJ.GT.MEM) GO TO 20
      IF(KL.LT.NROT) THEN
         KL=KL+1
         GO TO 10
      END IF
   20 CONTINUE
      MAXIJ=KL
      MEMUSE = IIJJ(MAXIJ)+MAXIJ
      IF(MEMUSE .GT. MEM) MAXIJ = MAXIJ-1
C
C     ----- SET TO ZERO -WXY- FOR THIS BLOCK -----
C
      DO 30 IJ=1,IIJJ(MAXIJ)+MAXIJ
   30    WXY(IJ) = ZERO
C
      CALL SEQREW(NFTI)
      CALL SEQREW(NFTA)
C
C     ----- (IJ|AB) INTEGRALS -----
C
      DO 60 J=1,NOC
      DO 60 I=J,NOC
         CALL TRFRD(NFTI,XX,NX,IX,NIX,NORB2)
         IF(DBUG) CALL TRFPRT(XX,IX,NORB,I,J)
         DO 60 K=1,NVIR
            IK=(K-1)*NOC+I
            IF(IK.GE.MINIJ.AND.IK.LE.MAXIJ) THEN
               KL=IA(K+NOC)+NOC
               DO 40 L=1,K
                  KL=KL+1
                  JL=(L-1)*NOC+J
   40             WXY(IIJJ(IK)+JL)=WXY(IIJJ(IK)+JL) - XX(KL)
            END IF
            JK=(K-1)*NOC+J
            IF(I.NE.J.AND.JK.GE.MINIJ.AND.JK.LE.MAXIJ) THEN
               KL=IA(K+NOC)+NOC
               DO 50 L=1,K-1
                  KL=KL+1
                  IL=(L-1)*NOC+I
   50             WXY(IIJJ(JK)+IL)=WXY(IIJJ(JK)+IL) - XX(KL)
            END IF
   60 CONTINUE
C
C     ----- (AI|BJ) INTEGRALS -----
C
      DO 100 J=1,NOC
      DO 100 I=1,NVIR
         JI=(I-1)*NOC+J
         CALL TRFRD(NFTA,XX,NX,IX,NIX,NORB2)
         IF(DBUG) CALL TRFPRT(XX,IX,NORB,I+NOC,J)
         IF(JI.GE.MINIJ.AND.JI.LE.MAXIJ) THEN
            DO 70 K=1,I-1
            JK=(K-1)*NOC+J
            KL=IA(K+NOC)
            DO 70 L=1,NOC
               KL=KL+1
               LK=(K-1)*NOC+L
   70          WXY(IIJJ(JI)+LK)=WXY(IIJJ(JI)+LK)+FOUR*XX(KL)
            KL=IA(I+NOC)
            DO 80 L=1,J
               KL=KL+1
               LI=(I-1)*NOC+L
   80          WXY(IIJJ(JI)+LI)=WXY(IIJJ(JI)+LI)+THREE*XX(KL)
         END IF
         DO 90 K=1,I-1
         KL=IA(K+NOC)
         JK=(K-1)*NOC+J
         DO 90 L=1,NOC
            LI=(I-1)*NOC+L
            KL=KL+1
            IF(LI.GE.MINIJ .AND. LI.LE.MAXIJ)
     *         WXY(IIJJ(LI)+JK)=WXY(IIJJ(LI)+JK)-XX(KL)
   90    CONTINUE
  100 CONTINUE
C
C     ----- ENERGY TERM TO DIAGONAL -----
C
      DO 110 I=1,NVIR
      DO 110 J=1,NOC
         JI=(I-1)*NOC+J
         IF(JI.GE.MINIJ.AND.JI.LE.MAXIJ)
     *      WXY(IIJJ(JI)+JI) = WXY(IIJJ(JI)+JI) + E(I+NOC) - E(J)
  110 CONTINUE
C
C     ----- WRITE TO FILE -----
C
      IF(NIRRED.EQ.1) THEN
         DO 120 JI=MINIJ,MAXIJ
            T(JI,5) = ONE/WXY(IIJJ(JI)+JI)
            CALL TRFWTM(NFT,WXY(IIJJ(JI)+1),NX,IX,NIX,JI,1,CUTOFF)
            NO = NO+NX
  120    CONTINUE
      ELSE
         I0=MINIJ/NOC+1
         J0=MOD(MINIJ,NOC)
         DO 140 I=I0,NVIR
         DO 140 J=J0,NOC
            J0=1
            IF(MOSIRP(I+NOC).NE.MOSIRP(J)) GO TO 140
            JI=(I-1)*NOC+J
            IF(JI.GT.MAXIJ) GO TO 150
            MROT = MROT+1
            T(MROT,5) = ONE/WXY(IIJJ(JI)+JI)
            KKLL=0
            KL=0
            DO 130 K=1,I
            L1=NOC
            IF(K.EQ.I) L1=J
            DO 130 L=1,L1
               KL=KL+1
               IF(MOSIRP(K+NOC).EQ.MOSIRP(L)) THEN
                  KKLL=KKLL+1
                  WXY(IIJJ(JI)+KKLL) = WXY(IIJJ(JI)+KL)
               END IF
  130       CONTINUE
            CALL TRFWTM(NFT,WXY(IIJJ(JI)+1),NX,IX,NIX,KKLL,1,CUTOFF)
            NO=NO+NX
  140    CONTINUE
  150 CONTINUE
      END IF
C                          DONE WITH THIS SLICE OF THE HESSIAN
      MINIJ=MAXIJ+1
      IF(MINIJ.LE.NROT) GO TO 1000
C
C     ----- SOLVE THE MPCPHF EQUATIONS -----
C
C     ----- FIRST COPY THE LAGRANGIAN TO -B- -----
C
      IF(NIRRED.EQ.1) THEN
         MROT=NROT
         IJ=0
         DO 200 I=1,NVIR
         DO 200 J=1,NOC
            IJ=IJ+1
  200       B(IJ) = XLAI(I,J)*T(IJ,5)
      ELSE
         MROT = 0
         IJ   = 0
         DO 210 I=1,NVIR
         DO 210 J=1,NOC
         IF(MOSIRP(I+NOC).EQ.MOSIRP(J)) THEN
            MROT=MROT+1
            B(MROT) = XLAI(I,J)*T(MROT,5)
         END IF
  210    CONTINUE
      END IF
C
C     ----- SOLVE THE SYSTEM  - WXY * XX = B - -----
C
      IF(SOME) WRITE(IW,9996) MROT,NO,NFT
      IF(MEM.GE.MROT+2*NO) THEN
         IF(SOME) WRITE(IW,9995) MROT+2*NO
         CALL MPCPRD(NFT,XX,IX,WXY,WXY(MROT+1),WXY(MROT+NO+1),
     *               MROT,CUTOFF)
         J01=J00+MROT
         J00=(J00-1)*NWDVAR+1
         J01=(J01-1)*NWDVAR+1
         DO 220 I=1,MROT
  220       XX(I) = ZERO
         CALL SPARSQ(WXY(MROT+NO+1),NROT,B,XX,MROT,T(1,1),T(1,2),
     *               T(1,3),T(1,4),MPSPAX,MPSPAX)
      ELSE
         DO 235 I=1,MROT
  235       XX(I) = ZERO
         CALL SPARSQ(WXY,NROT,B,XX,MROT,T(1,1),T(1,2),T(1,3),T(1,4),
     *               MPOUAX,MPOUAX)
      END IF
C
C     ----- COPY THE SOLUTION BACK TO -XLAI- -----
C
      IF(NIRRED.EQ.1) THEN
         IJ = 0
         DO 240 I=1,NVIR
         DO 240 J=1,NOC
            IJ=IJ+1
  240       XLAI(I,J) = XX(IJ)
      ELSE
         IJ=0
         KL=0
         DO 250 I=1,NVIR
         DO 250 J=1,NOC
            IJ=IJ+1
            IF(MOSIRP(I+NOC).EQ.MOSIRP(J)) THEN
               KL=KL+1
               XLAI(I,J) = XX(KL)
            ELSE
               XLAI(I,J) = ZERO
            END IF
  250    CONTINUE
      END IF
C
      IF(DBUG) THEN
         WRITE(IW,9997)
         CALL PRSQ(XLAI,NOC,NVIR,NVIR)
      END IF
C
      RETURN
C
 9997 FORMAT(/' PAI(VIRT,OCC)')
 9996 FORMAT(5X,'CHOOSING OUT OF CORE SOLUTION OF THE CPHF EQUATIONS,'/
     *       5X,'THE SIZE OF THE LINEAR EQUATION SYSTEM IS',I10/
     *       5X,I8,' ELEMENTS OF THE ORBITAL HESSIAN ARE STORED ON',
     *             ' FILE',I3)
 9995 FORMAT(5X,'USING FAST PACKED ROUTINES WITH NWORD=',I10)
      END
C
C*MODULE MP2GRD  *DECK MPCPRD
      SUBROUTINE MPCPRD(NFTI,XX,IX,II,JJ,WXY,NROT,TOL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION XX(NROT),IX(NROT),II(NROT),JJ(*),WXY(*)
C
      CALL SEQREW(NFTI)
      IJ = 0
      DO 20 IROT=1,NROT
         CALL TRFRD(NFTI,XX,NX,IX,NIX,NROT)
         ICOUNT = 0
         DO 10 JROT=1,IROT
            IF(ABS(XX(JROT)).GT.TOL) THEN
               IJ=IJ+1
               WXY(IJ) = XX(JROT)
               JJ(IJ) = JROT
               ICOUNT = ICOUNT+1
            END IF
   10    CONTINUE
         II(IROT) = ICOUNT
   20 CONTINUE
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MPGTRP
      SUBROUTINE MPGTRP(PIJ,PAB,PAI,C,NBF,NOC,NVIR,
     *           X,Y,PMN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00,TWO=2.0D+00)
C
      DIMENSION PIJ(NOC,NOC),PAB(NVIR,NVIR),PAI(NVIR,NOC),
     *          C(NBF,*),X(NBF,*),Y(NBF,*),PMN(*)
C
      CALL MRARTR(PIJ,NOC,NOC,NOC,C,NBF,NBF,X,NBF)
      CALL MRARBR(C,NBF,NBF,NOC,X,NBF,NBF,Y,NBF)
C
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,ONE,
     *     Y,NBF)
C
      CALL MRARTR(PAI,NVIR,NVIR,NOC,C,NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,TWO,C(1,NOC+1),NBF,X,NBF,ONE,
     *     Y,NBF)
C
      MUNU=0
      DO 15 MU=1,NBF
         DO 10 NU=1,MU
            MUNU=MUNU+1
            PMN(MUNU) = Y(MU,NU) + Y(NU,MU)
   10    CONTINUE
   15 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MPOUAX
      SUBROUTINE MPOUAX(A,LDA,X,AX,N)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /MP2PTR/ NFT,J00,J01,J02,IDUMMY(4)
      COMMON /FMCOM / XX(1)
C
      DIMENSION A(*),X(*),AX(*)
C
C     ----- THIS ROUTINE FORM THE PRODUCT OF A MATRIX -A- AND
C           A VECTOR -X- WITH THE MATRIX -A- STORED IN DISK
C           THE MATRIX IS SUPPOSED TO BE SYMMETRIC
C
      CALL SEQREW(NFT)
      DO 10 I=1,N
   10    AX(I)=ZERO
C
      DO 40 I=1,N
         CALL TRFRD(NFT,A,NX,XX(J01),NIX,N)
         DO 20 J=1,I
   20       AX(I)=AX(I) + A(J)*X(J)
         DO 30 J=1,I-1
   30       AX(J)=AX(J) + A(J)*X(I)
   40 CONTINUE
C
      DO 50 I=1,N
   50    AX(I)=AX(I)*XX(J02+I-1)
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MPSPAX
      SUBROUTINE MPSPAX(A,LDA,X,AX,N)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /MP2PTR/ NFT,J00,J01,J02,IDUMMY(4)
      COMMON /FMCOM / XX(1)
C
      DIMENSION IX(1)
      EQUIVALENCE (XX(1),IX(1))
      DIMENSION A(*),X(*),AX(*)
C
      II(I) = IX(J00-1+I)
      JJ(I) = IX(J01-1+I)
C
      DO 10 I=1,N
         AX(I) = ZERO
   10 CONTINUE
C
      IJ = 0
      DO 40 I=1,N
         NUM = II(I)
         DO 20 J=1,NUM
            AX(I) = AX(I) + A(IJ+J)*X(JJ(IJ+J))
   20    CONTINUE
         DO 30 J=1,NUM-1
            AX(JJ(IJ+J)) = AX(JJ(IJ+J)) + A(IJ+J)*X(I)
   30    CONTINUE
         IJ = IJ+NUM
   40 CONTINUE
C
      DO 50 I=1,N
         AX(I)=AX(I)*XX(J02-1+I)
   50 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MPW2IJ
      SUBROUTINE MPW2IJ(POCC,PVIR,PAI,WOCC,WAI,E,XX,
     *                  X,NOC,NVIR,NORB,NORB2,NFTI,NFTA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
      DIMENSION POCC(NOC,NOC),PVIR(NVIR,NVIR),PAI(NVIR,NOC),
     *          WOCC(NOC,NOC),WAI(NVIR,NOC)
      DIMENSION X(NORB,NORB),XX(NORB2),E(*)
C
      CALL SEQREW(NFTI)
      CALL SEQREW(NFTA)
C
      DO 75 J=1,NOC
         DO 70 I=J,NOC
            CALL TRFRD(NFTI,XX,NX,X,NIX,NORB2)
            KL=0
            DO 15 K=1,NORB
               DO 10 L=1,K
                  KL=KL+1
                  X(K,L) = XX(KL)
                  X(L,K) = XX(KL)
   10          CONTINUE
   15       CONTINUE
C
C ... (IJ|KL) INTEGRALS
C
         TMP1 = ZERO
         DO 25 L=1,NOC
            DO 20 K=1,NOC
               TMP1 = TMP1 + POCC(K,L)*X(K,L)
   20       CONTINUE
   25    CONTINUE
C
         DO 35 L=1,NOC
            DO 30 K=1,NOC
               WOCC(I,L) = WOCC(I,L) + POCC(K,J)*X(K,L)
   30       CONTINUE
   35    CONTINUE
C
         IF(I.NE.J) THEN
            DO 45 L=1,NOC
               DO 40 K=1,NOC
                  WOCC(J,L) = WOCC(J,L) + POCC(K,I)*X(K,L)
   40          CONTINUE
   45       CONTINUE
         END IF
C
C ... (IJ|AL) INTEGRALS
C
         TMP2 = ZERO
         DO 55 L=1,NOC
            DO 50 K=1,NVIR
               TMP2 = TMP2 + PAI(K,L)*X(K+NOC,L)
   50       CONTINUE
   55    CONTINUE
         TMP1 = TMP1 + TWO*TMP2
C
C ... (IJ|AB) INTEGRALS
C
         DO 65 L=1,NVIR
            DO 60 K=1,NVIR
               TMP1 = TMP1 + PVIR(K,L)*X(K+NOC,L+NOC)
   60       CONTINUE
   65    CONTINUE
C
         WOCC(I,J) = WOCC(I,J) - TWO*TMP1
         IF(I.NE.J) WOCC(J,I) = WOCC(J,I) - TWO*TMP1
C
   70    CONTINUE
   75 CONTINUE
C
      DO 115 J=1,NOC
         DO 110 I=1,NVIR
            CALL TRFRD(NFTA,XX,NX, X,NIX,NORB2)
            KL=0
            DO 85 K=1,NORB
               DO 80 L=1,K
                  KL=KL+1
                  X(K,L) = XX(KL)
                  X(L,K) = XX(KL)
   80          CONTINUE
   85       CONTINUE
C
C ... (AI|KL) INTEGRALS
C
            DO 95 L=1,NOC
               DO 90 K=1,NOC
                  WOCC(L,J) = WOCC(L,J) + TWO*PAI(I,K)*X(K,L)
   90          CONTINUE
   95       CONTINUE
C
C ... (AI|BJ) INTEGRALS
C
            DO 105 L=1,NOC
               DO 100 K=1,NVIR
                  WOCC(L,J) = WOCC(L,J) + PVIR(K,I)*X(K+NOC,L)
  100          CONTINUE
  105       CONTINUE
C
  110    CONTINUE
  115 CONTINUE
C
C ... ENERGY TERM TO WAI
C
      DO 125 I=1,NOC
         DO 120 J=1,NVIR
            WAI(J,I) = WAI(J,I) - E(I)*PAI(J,I)
  120    CONTINUE
  125 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MPWMP2
      SUBROUTINE MPWMP2(V,E,D1,D2,NOC,N1,KREC7,KREC8,KREC9)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C     COMMON /DAFILE/ IDAF,NAV,IODA(255)
C
      DIMENSION D1(*),D2(*),V(N1,*),E(*)
C
      N2=(N1*N1+N1)/2
C
C     -----  D2 IS THE HF+MP2 ENERGY-WEIGHTED DENSITY  -----
C
      OCC=TWO
      IJ=0
      DO 25 I=1, N1
         DO 20 J=1, I
            IJ=IJ+1
            DUM=ZERO
            DO 30 K=1, NOC
               DUM=DUM-E(K)*V(I,K)*V(J,K)
   30       CONTINUE
            D2(IJ)=D2(IJ)+DUM*OCC
   20    CONTINUE
   25 CONTINUE
C
      CALL DAWRIT(IDAF,IODA,D2,N2,KREC9,0)
C
C     -----  FORM TOTAL HF+MP2 DENSITY  -----
C
      CALL DAREAD(IDAF,IODA,D1,N2,16,0)
      CALL DAREAD(IDAF,IODA,D2,N2,KREC7,0)
      DO 10 N=1, N2
         D2(N)=D1(N) + D2(N)
   10 CONTINUE
C
      CALL DAWRIT(IDAF,IODA,D2,N2,16,0)
C
C     ----- SAVE PHF ON KREC8 -----
C
      CALL DAWRIT(IDAF,IODA,D1,N2,KREC8,0)
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MP2GC2
      SUBROUTINE MP2GC2(XPQRJ,C,X,Y,XX,MOSIRP,NOC,NBF,NBF2,NORB,TOL,
     *                  NO2EA,NO2EI,NOCC1,NOCC,NFTIA,NFTIJ,NCO,MOCPHF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      LOGICAL SOME,OUT,DBUG,SYM2EI
C
      COMMON /MP2PRT/ SOME,OUT,DBUG
C
      DIMENSION XPQRJ(NOC,NBF2,NBF),C(NBF,NORB),
     *          X(NBF,*),Y(NBF,*),XX(*),MOSIRP(NORB)
C
      NORB2 = (NORB*NORB+NORB)/2
C
C     FINAL THREE INDICES OF THE INTEGRAL TRANSFORMATION
C     NOTE: VARIABLES HAVE CHANGED SINCE MP2GE2
C
C     ----- (PQ|RJ) TO (PQ|AJ) -----
C     IF WE ARE DOING CPHF EQUATIONS IN AO BASIS, WE ONLY NEED TO DO
C     THE SECOND TRANSFORMATION OVER THE VIRTUAL ORBITALS.  IF WE WILL
C     DO THE CPHF EQUATIONS IN THE MO BASIS, WE NEED TO DO THE SECOND
C     TRANSFORMATION OVER THE ENTIRE SET OF ORBITALS AND SAVE THE
C     (IJ|**) SET ON DISK.
C
      IF (MOCPHF.EQ.0) THEN
         NP1=NCO+NOCC+1
         NTRAN=NORB-NCO-NOCC
      ELSE
         NP1=1
         NTRAN=NORB
      END IF
C
      DO 20 MPQ=1,NBF2
         DO 15 MR=1,NBF
            DO 10 MJ=1,NOC
               X(MR,MJ) = XPQRJ(MJ,MPQ,MR)
   10       CONTINUE
   15    CONTINUE
C
         CALL MRTRBR(C(1,NP1),NBF,NBF,NTRAN,X,NBF,NOC,Y(NP1,1),NBF)
C
         DO 30 MA=1,NTRAN
            DO 25 MJ=1,NOC
               XPQRJ(MJ,MPQ,MA+NP1-1) = Y(MA+NP1-1,MJ)
   25       CONTINUE
   30    CONTINUE
   20 CONTINUE
C
C     ----- (PQ|AJ) TO (BC|AJ) -----
C
      DO 80 MJ=1,NOC
      DO 75 MA=1,NORB
         IF(MA.LT.MJ+NOCC1-1+NCO) GO TO 75
         JKT = NFTIA
         IF(MA.LE.NOCC+NCO) THEN
            IF(MOCPHF.EQ.0) GO TO 75
            IF(MOCPHF.EQ.1) JKT = NFTIJ
         END IF
         MPQ = 0
         DO 40 MP=1,NBF
            DO 35 MQ=1,MP
               MPQ = MPQ+1
               X(MP,MQ) = XPQRJ(MJ,MPQ,MA)
               X(MQ,MP) = X(MP,MQ)
   35       CONTINUE
   40    CONTINUE
         CALL MRARBR(X,NBF,NBF,NBF,C,NBF,NORB,Y,NBF)
         MPQ=0
         DO 50 MP=1,NORB
            DO 45 MQ=1,MP
               MPQ = MPQ+1
               IF(SYM2EI(MOSIRP(MP),MOSIRP(MQ),
     *                   MOSIRP(MA),MOSIRP(MJ+NCO+NOCC1-1))) THEN
                  XX(MPQ) = DDOT(NBF,Y(1,MP),1,C(1,MQ),1)
                  IF(ABS(XX(MPQ)).LT.TOL) XX(MPQ) = ZERO
               ELSE
                  XX(MPQ) = ZERO
               END IF
   45       CONTINUE
   50    CONTINUE
C
         DO 60 MPQ=1,NORB2
            XPQRJ(MJ,MPQ,MA) = XX(MPQ)
   60    CONTINUE
C
         IF(MA.GT.NOCC+NCO) THEN
            MPQ=0
            DO 70 MP=1,NORB
               DO 65 MQ=1,MP
                  MPQ=MPQ+1
                  IF(MP.GT.NOCC+NCO.AND.MQ.GT.NOCC+NCO) XX(MPQ) = ZERO
   65          CONTINUE
   70       CONTINUE
         END IF
C
         CALL TRFWTM(JKT,XX,NX,Y,NIX,NORB2,1,TOL)
         IF(DBUG.AND.NX.GT.0) THEN
            CALL TRFOUT(XX,NX,Y,NIX,NORB,MA,MJ+NOCC1-1+NCO)
         END IF
         IF(JKT.EQ.NFTIA) THEN
            NO2EA = NO2EA+NX
         ELSE
            NO2EI = NO2EI+NX
         END IF
C
   75 CONTINUE
   80 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MP2GE2
      SUBROUTINE MP2GE2
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXRT=100, MXATM=2000, MXAO=8192)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,DIRSCF,FDIFF,GOPARR,DSKWRK,MASWRK,
     *        ABEL,ABELPT,DIRTRF
C
      COMMON /DIRMEM/ IPQKL,IVEC,IDUM1,IDUM2,NUM,NOCC,NDIM,MINPQ,MAXPQ
      COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /ENRGMP/ EMP2,EMP3,EMP4
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NNUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWORD,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /MP2PRT/ SOME,OUT,DBUG
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      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 /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
C     COSMO ADDITIONS FOR MP2 ANALYSIS
C
      LOGICAL ISEPS, USEPS
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /COSMO2/ QENUC,ELAST,EMP2COS,EMP2COS2,ETOTSMP,SAVESE,
     *                EMP2LAST,MP2TRIP,MP2ITER,MP2FACT
C
      DATA KREC1,KREC2,KREC3,KREC4,KREC5,KREC6/301,302,303,304,305,306/
      DATA NFTIA,NFTIJ,NFT16/9,15,16/
      DATA CHECK/8HCHECK   /, DBUGME,DEBUG/8HMP2GRAD ,8HDEBUG   /
C
C     ----- THIS ROUTINE CALCULATES THE ENERGY CONTRIBUTIONS TO THE
C     MP2 GRADIENT.  AFTER THE ENERGY IS CALCULATED, SOME PART OF
C     DENSITY AND ENERGY WEIGHTED DENSITY ARE CALCULATED -----
C
C     ----- THIS CODE IS WRITTEN FROM THE PAPER BY M. DUPUIS, S. CHIN
C     AND A. MARQUEZ IN RELATIVISTIC AND ELECTRON CORRELATION EFFECTS
C     IN MOLECULES AND SOLIDS, G. L. MALI ED., PLENUM PRESS, NEW YORK,
C     1994, P 315.  AND THE ABOVE IS BASED ON THE PAPER M. J. FRISCH,
C     M. HEAD-GORDON, J. A. POPLE CHEM. PHYS. LETT.166 (1990) 275 -----
C
C     ---- MODIFICATIONS TO ENABLE USE OF THE FROZEN CORE APPROXIMATION
C     WERE MADE BY S.P.WEBB AND WERE COMPLETED JUNE 1996 AT ISU ----
C
      CUTOFF=TOL
      NOCA  = NOA - NACORE
      NVIRA = NO  - NOA
      NORB  = NO  - NACORE
      NOC   = NOA
      NVIR  = NVIRA
C
      NBF2  = (NBF*NBF+NBF)/2
      NBF3  = NBF*NBF
      NUM   = NBF
      NDIM  = NBF2
C
      SOME = NPRINT.NE.-5  .AND.  MASWRK
      OUT  = .FALSE.
      DBUG = MASWRK  .AND.  (EXETYP.EQ.DBUGME  .OR.  EXETYP.EQ.DEBUG)
      IF(MASWRK) WRITE(IW,9000)
      IF(SOME) WRITE(IW,9005) NACORE,NOCA,NQMT,NBF
C
      CALL SEQOPN(NFTIA,'MOINTS','UNKNOWN',.FALSE.,'UNFORMATTED')
      IF (MOCPHF.EQ.1)
     *  CALL SEQOPN(NFTIJ,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFT16,'WORK16','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFTIA)
      IF (MOCPHF.EQ.1) CALL SEQREW(NFTIJ)
C
C        ---- OPTION TO DEFEAT SYMMETRY USAGE ----
C
      ABEL=ABELPT()
      IF(NOSYMT.EQ.1) CALL SYMOFF
C
      CALL GOTFM(NGOTMX)
      IF(NWORD.GT.0) NGOTMX=MIN(NWORD,NGOTMX)
C
      CALL VALFM(LOADFM)
      LVEC  = 1     + LOADFM
      IIRP  = LVEC  + NBF3
      LAST  = IIRP  + NBF
      LVEC2 = LVEC
      IF(NACORE.GT.0) THEN
        LVEC2 = LAST
        LAST  = LVEC2 + NBF3
      END IF
      NEEDA = LAST  - LOADFM
      CALL GETFM(NEEDA)
C
C     ----- READ IN ORBITALS, AND ASSIGN SYMMETRY -----
C
      CALL VALFM(LOADFM)
      IWRK   = 1     + LOADFM
      IWRK2  = IWRK  + NBF
      ITVEC  = IWRK2 + NBF
      IPOCC  = ITVEC + NBF2
      ITMP   = IPOCC + NBF3
      LAST   = ITMP  + NBF
      NEEDB  = LAST  - LOADFM
      CALL GETFM(NEEDB)
      IF(EXETYP.EQ.CHECK) GO TO 90
C
      IF (NACORE.GT.0) THEN
C
C     ---- NEED ONLY ACTIVE ORBITAL COEFFICIENTS IN 1 INDEX ----
C     ---- TRANSFORMATION (INPENDING CALL TO MP2C1)         ----
C
         CALL DAREAD(IDAF,IODA,X(IPOCC) ,NBF3,15,0)
         DO 50 J=1,NORB
            DO 60 I=1,NBF
               X(I+ NBF*(J-1)+LVEC-1)=X(I+NBF*(J-1+NACORE)+IPOCC-1)
  60        CONTINUE
  50     CONTINUE
C
C     ---- NEED ALL ORBITAL COEFFICIENTS TO ASSIGN SYMMETRY ----
C     ---- AND TRANSFORMATION OF REMAINING INDICES (MP2GC2) ----
C
         CALL DAREAD(IDAF,IODA,X(LVEC2) ,NBF3,15,0)
C
         CALL DAREAD(IDAF,IODA,X(ITVEC),NBF2,12,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
C
         CALL TRFSYM(X(IWRK),X(IIRP),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVEC2),X(ITMP),IA,NQMT,NBF,NO,NBF)
C
      ELSE
C
         CALL DAREAD(IDAF,IODA,X(LVEC) ,NBF3,15,0)
         CALL DAREAD(IDAF,IODA,X(ITVEC),NBF2,12,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
         CALL TRFSYM(X(IWRK),X(IIRP),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVEC),X(ITMP),IA,NQMT,NBF,NO,NBF)
      END IF
C
   90 CONTINUE
      CALL RETFM(NEEDB)
C
C     -LVEC - -I00-  NBF*NBF           - ORBITAL COEFF
C     -IIRP - -I01-  NBF               - MOSIRP
C     -IWRK - -I02-  MAX(INITMX,NBF3)  - INTERMEDIATE DATA
C     -IWRK2- -I03-  MAX(NINTMX,NBF3)  - INTERMEDIATE DATA
C     -ITVEC- -I04-  NOC*NBF           - TRANSPOSE C
C     -IPOCC- -I05-  NOC*NOC           - OCC-OCC BLOCK DENSITY (POCC)
C     -IWOCC- -I06-  NOC*NOC           - OCC-OCC E WEIGHTED DENS (WOCC)
C     -IPVIR- -I07-  NVIR*NVIR         - VIR-VIR DENSITY (PVIR)
C     -IWVIR- -I08-  NVIR*NVIR         - VIR-VIR E WEIGHTED DENS (WVIR)
C     -ILAI1- -I09-  NOC*NVIR          - MP2 LAGRANGIAN LAI(I)
C     -ILAI2- -I10-  NOC*NVIR          - MP2 LAGRANGIAN LAI(II)
C     -IPINT- -I11-  NBF*NBF2*NOC      - (PQ|RJ) INTEGRALS
C     -IWRK3- -I12-  NOC*NVIR*NOC      - INTERMEDIATE DATA
C
      CALL VALFM(LOADFM)
      IWRK   = 1     + LOADFM
      IWRK2  = IWRK  + MAX(NINTMX,NBF3)
      ITVEC  = IWRK2 + MAX(NINTMX,NBF3)
      IPOCC  = ITVEC + NBF
      IWOCC  = IPOCC + NOC*NOC
      IPVIR  = IWOCC + NOC*NOC
      IWVIR  = IPVIR + NVIR*NVIR
      ILAI1  = IWVIR + NVIR*NVIR
      ILAI2  = ILAI1 + NOC*NVIR
      LAST   = ILAI2 + NOC*NVIR
      LSTSAV = LAST
C
C        FOR THIS I2CASE=3 DIRECT TEI TRANSFORMATION, THE /DIRMEM/
C        COMMON MUST INITIALIZE -IPQKL-, -IVEC- (POINTING AT FIRST
C        MO OF EACH PARTICULAR PASS), -NUM-, AND -NOCC- (WHICH IS
C        SET INSIDE THE MP2C1 ROUTINE).  OTHER ELEMENTS ARE NOT USED.
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
      NSH2 = (NSHELL*NSHELL+NSHELL)/2
      IF(DIRSCF) THEN
         IGHND  = LAST   + NOC*NVIR
         IXINTS = IGHND  + MAXG
         IDDIJ  = IXINTS + NSH2
         LAST   = IDDIJ  + 49*MXG2
      END IF
      NEEDC  = LAST - LOADFM
      CALL GETFM(NEEDC)
C
      IF(DIRSCF  .AND.  ISCHWZ.EQ.1)
     *     CALL DAREAD(IDAF,IODA,X(IXINTS),NSH2,54,0)
C
      LEFT   = NGOTMX - NEEDA - NEEDC
      NMIN   = NEEDA + NEEDC + NVIR*NOCA     + NBF2*NBF
      NMAX   = NEEDA + NEEDC + NVIR*NOCA*NOCA + NBF2*NBF*NOCA
      NINCR  =                 NVIR*NOCA    + NBF2*NBF
      IF(SOME) WRITE(IW,9990) NMIN,NMAX,NINCR
C
      IF(NMIN.GT.NGOTMX) THEN
         IF(MASWRK) WRITE(IW,9997) NMIN-LEFT
         CALL ABRT
      END IF
C
      LPASS = MIN(NOCA,LEFT/NINCR)
      NPASS = (NOCA+LPASS-1)/LPASS
      LPASS = (NOCA+NPASS-1)/NPASS
C
      NOCC1  = 1
      NO2EI  = 0
      NO2EA  = 0
C
      DO 102 I=IPOCC,LSTSAV
  102 X(I)=ZERO
C
      CALL VALFM(LOADFM)
      IPQKL  = 1     + LOADFM
      IWRK3  = IPQKL + NBF*NBF2*LPASS
      LAST   = IWRK3 + NVIR*NOC*LPASS
      NEEDD  = LAST  - LOADFM
      CALL GETFM(NEEDD)
      NTOT  = NEEDA + NEEDC + NEEDD
      IF(SOME) WRITE(IW,9999) NGOTMX,NTOT,LPASS,NPASS,CUTOFF
      IF(EXETYP.EQ.CHECK) GO TO 600
C
      IPASS = 0
      CALL TSECND(TIM0)
  100 CONTINUE
      IPASS = IPASS+1
      NOCC2  = MIN(NOCA,NOCC1+LPASS-1)
      NOCC   = NOCC2-NOCC1+1
      IVEC   = LVEC+(NOCC1-1)*NBF
C
C     ----- TRANSFORMATION, THEN FILL PART OF DENSITY MATRICES -----
C
      NTSAVE = NT
      IF(DIRSCF  .AND.  .NOT.ABEL) NT=1
C
      CALL MP2C1(IJK,X(IPQKL),X(IWRK),X(IWRK2),X(IVEC),
     *           NBF,NOCC,NBF2,X(IGHND),X(IXINTS),X(IDDIJ),MAXG,NSH2)
      CALL MP2GC2(X(IPQKL),X(LVEC2),X(IWRK),X(IWRK2),X(IWRK),
     *            X(IIRP),NOCC,NBF,NBF2,NO,CUTOFF,NO2EA,NO2EI,
     *            NOCC1,NOCA,NFTIA,NFTIJ,NACORE,MOCPHF)
C
      NT = NTSAVE
      IF(NOSYMT.EQ.1) CALL SYMON
C
      CALL DAREAD(IDAF,IODA,X(ITVEC),NBF,17,0)
      CALL P2W2N(X(IPQKL),X(IWRK3),CUTOFF,X(ITVEC),X(IPOCC),
     *           X(IWOCC),X(IPVIR),X(IWVIR),X(ILAI1),X(ILAI2),
     *           NORB,NOC,NVIR,NOCC,NOCC1,NBF2,NBF,NACORE,NOCA)
C
      NOCC1 = NOCC2+1
      CALL TSECND(TIM1)
      TPASS = TIM1-TIM0
      TIM0  = TIM1
      IF(SOME) WRITE(IW,9980) IPASS,TPASS
      IF(SOME) CALL FLSHBF(IW)
      IF(NOCC1.LE.NOCA) GO TO 100
C
C     ----- GET MP2 ENERGY AND ADD ENERGY CONTRIBUTIONS TO
C           -WOCC- AND -WVIR- -----
C
      CALL MP2GEN(E2,X(ITVEC),X(IPOCC),X(IWOCC),X(IPVIR),X(IWVIR),
     *            NOC,NVIR)
C
      EHF  = ESCF
      EMP2 = EHF+E2
C
C     COSMO CHANGES
C
      IF(ISEPS) THEN
         WRITE(IW,*)"INSIDE MP2GE2:"
         WRITE(IW,*)"VALUE OF ISEPS IN MP2:",ISEPS
         IF(MP2TRIP.EQ.0) EMP2COS=EMP2
         MP2TRIP=1
         WRITE(IW,7802) EHF,E2,EMP2
      ELSE
         WRITE(IW,7800) EHF,E2,EMP2
      ENDIF
C
C
      IF(DBUG) THEN
         WRITE(IW,8999)
         CALL PRSQ(X(IPOCC),NOC,NOC,NOC)
         WRITE(IW,8998)
         CALL PRSQ(X(IPVIR),NVIR,NVIR,NVIR)
         WRITE(IW,8997)
         CALL PRSQ(X(IWOCC),NOC,NOC,NOC)
         WRITE(IW,8996)
         CALL PRSQ(X(IWVIR),NVIR,NVIR,NVIR)
         WRITE(IW,8995)
         CALL PRSQ(X(ILAI1),NOC,NVIR,NVIR)
         WRITE(IW,8994)
         CALL PRSQ(X(ILAI2),NOC,NVIR,NVIR)
      END IF
C
      IF(SOME) WRITE(IW,9996) NO2EI,NFTIJ,NO2EA,NFTIA
C
C     ----- SAVE ON DIRECT ACCESS FILE -----
C
C     -POCC- ON RECORD -KREC1-
C     -PVIR- ON RECORD -KREC2-
C     -WAI - ON RECORD -KREC3-
C     -LAI - ON RECORD -KREC4-
C     -WOCC- ON RECORD -KREC5-
C     -WVIR- ON RECORD -KREC6-
C         IN OTHER PLACES, WE WILL SAVE THREE OTHER RECORDS:
C     -P(2)- ON RECORD -KREC7-
C     -PHF-  ON RECORD -KREC8-
C     -W(2)- ON RECORD -KREC9-
C
  600 CONTINUE
      CALL DAWRIT(IDAF,IODA,X(IPOCC),NOC *NOC ,KREC1,0)
      CALL DAWRIT(IDAF,IODA,X(IPVIR),NVIR*NVIR,KREC2,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI2),NVIR*NOC ,KREC3,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI1),NVIR*NOC ,KREC4,0)
      CALL DAWRIT(IDAF,IODA,X(IWOCC),NOC *NOC ,KREC5,0)
      CALL DAWRIT(IDAF,IODA,X(IWVIR),NVIR*NVIR,KREC6,0)
C
C     ----- ENERGY IS CALCULATED, RETURN TO THE CALLING PROGRAM -----
C
      CALL RETFM(NEEDD)
      CALL RETFM(NEEDC)
      CALL RETFM(NEEDA)
      RETURN
C
 8999 FORMAT(/' ----- P2(OCC,OCC)     -----'/)
 8998 FORMAT(/' ----- P2(VIR,VIR)     -----'/)
 8997 FORMAT(/' ----- W2(OCC,OCC)     -----'/)
 8996 FORMAT(/' ----- W2(VIR,VIR)     -----'/)
 8995 FORMAT(/' ----- LAGRANGIAN (I)  -----'/)
 8994 FORMAT(/' ----- W2(VIR,OCC) (I) -----'/)
 9000 FORMAT(/1X,16('-'),3X,52(1H-)/
     *        1X,'RHF-MP2 GRADIENT',3X,
     *           'PROGRAM WRITTEN BY A.MARQUEZ, S.MALUENDES, M.DUPUIS,'/
     *       20X,'AND S.P.WEBB'/
     *        1X,16('-'),3X,52(1H-))
 9005 FORMAT(/1X,'# CORE ORBITALS      =',I5/
     *        1X,'# OCCUPIED ORBITALS  =',I5/
     *        1X,'# MOLECULAR ORBITALS =',I5/
     *        1X,'# BASIS FUNCTIONS    =',I5)
 7800 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/)
 7802 FORMAT(/1X,'RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE'/
     *       13X,'  ECOSMO(0)               =',1X,F20.10/
     *       13X,'  E(1)                    =',9X,'0.0'/
     *       13X,'  E(2)                    =',1X,F20.10/
     *       13X,'E(COSMO-MP2)              =',1X,F20.10/)
C
 9999 FORMAT(/' MP2 ENERGY/GRADIENT TEI TRANSFORMATION'/
     *        ' # OF WORDS AVAILABLE      = ',I10/
     *        ' # OF WORDS USED           = ',I10/
     *        ' # OF ORBITALS/PASS        = ',I10/
     *        ' # OF INTEGRAL PASSES      = ',I10/
     *        ' THRESHOLD FOR KEEPING TRANSFORMED 2E- INTEGRALS = ',
     *        1P,E10.3/)
 9997 FORMAT(/' NOT ENOUGH MEMORY IN -MP2- TRANSFORMATION'/
     *        ' INCREASE MEMORY SIZE BY AT LEAST ',I10,' WORDS')
 9996 FORMAT(1X,'DONE WITH -MP2- GRADIENT TEI TRANFORMATION'/
     *       1X,I10,' NON-ZERO [IJ|AB] INTEGRALS WRITTEN TO FILE',I3/
     *       1X,I10,' NON-ZERO [AI|BJ] INTEGRALS WRITTEN TO FILE',I3)
 9980 FORMAT(1X,'PASS #',I4,' TOOK',F10.2,' SECONDS.')
 9990 FORMAT(/1X,'MEMORY REQUIREMENTS FOR MP2 TRANSFORMATION'/
     *       1X,' MINIMUM=',I10,' WORDS, USING 1 ORBITAL PER PASS'/
     *       1X,' MAXIMUM=',I10,' WORDS, USING ONLY 1 AO INTEGRAL PASS'/
     *       1X,'EACH ADDITIONAL',I10,' WORDS REDUCES PASS COUNT BY 1.')
      END
C
C*MODULE MP2GRD  *DECK MP2GEN
      SUBROUTINE MP2GEN(EMP2,E,POCC,WOCC,PVIR,WVIR,NOC,NVIR)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00)
C
      DIMENSION E(*),POCC(NOC,NOC),WOCC(NOC,NOC),PVIR(NVIR,NVIR),
     *          WVIR(NVIR,NVIR)
C
C     ----- MP2 ENERGY -----
C
      EMP2 = ZERO
      DO 10 MI=1,NOC
         EMP2 = EMP2 + WOCC(MI,MI)
   10 CONTINUE
      EMP2 = -EMP2
C
C     ----- ENERGY CONTRIBUTIONS TO -WOCC- AND -WVIR- -----
C
      DO 25 MJ=1,NOC
         DO 20 MI=1,NOC
            WOCC(MI,MJ) = WOCC(MI,MJ) - POCC(MI,MJ)*(E(MI)+E(MJ))*PT5
   20    CONTINUE
   25 CONTINUE
C
      DO 35 MB=1,NVIR
         DO 30 MA=1,NVIR
            WVIR(MA,MB) = WVIR(MA,MB) -
     *                    PVIR(MA,MB)*(E(MA+NOC)+E(MB+NOC))*PT5
   30    CONTINUE
   35 CONTINUE
      RETURN
      END
C
C*MODULE MP2GRD  *DECK MP2GRD
      SUBROUTINE MP2GRD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000, MXRT=100, MXAO=8192)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,GOPARR,DSKWRK,MASWRK,
     *        ABEL,ABELPT,DIRSCF,FDIFF,DIRTRF
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      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 /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
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWORD,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /MP2PTR/ NFT,J00,J01,J02,IDUMMY(4)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFMEM/ IFCK,IDUM,IDEN,JDUM(3),IBUF,KDUM(3)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
C     COSMO INFORMATION
C
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
      LOGICAL ISEPS, USEPS
      DIMENSION QVPOT(NPPA)
      COMMON /ISEPS/  ISEPS, USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     1                DEN(MAXDEN),AR(LENABC),
     3                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /SOLVI / IATSP(LENABC+1),N0(2),NP1,NP2,ISKIP
      COMMON /COSMO1/ SE2,SECORR,ETOTS,CDUM,QVCOSMO,
     *                CSPOT(NPPA),ICORR,ITRIPO,ITRIP2,ITRIP3,ITRIP4,
     *                NATCOS,NQS,ITERC
      COMMON /COSMO2/ QENUC,ELAST,EMP2COS,EMP2COS2,ETOTSMP,SAVESE,
     *                EMP2LAST,MP2TRIP,MP2ITER,MP2FACT
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      DATA KREC1,KREC2,KREC3,KREC4,KREC5,KREC6/301,302,303,304,305,306/
      DATA KREC7,KREC8,KREC9/307,308,309/
      DATA NFTIA,NFTIJ,NFT16/9,15,16/
      DATA CHECK/8HCHECK   /
C
C     ----- THIS ROUTINE COMPLETES THE SETUP FOR THE MP2 GRADIENT -----
C     FIRST, FORM MP2 LAGRANGIAN, AND SOLVE THE CPHF-LIKE EQUATION.
C     THEN, FORM THE FIRST ORDER WAVEFUNCTION'S FIRST ORDER DENSITY,
C     AND THE SO-CALLED ENERGY-WEIGHTED DENSITY MATRIX.
C     AT THIS POINT, IT IS POSSIBLE TO DO 1E- PROPERTIES.
C     FINALLY, CREATE THE NON-SEPARABLE PART OF THE 2ND ORDER DENSITY.
C
      SOME = MASWRK  .AND.  NPRINT.NE.-5
      IF(SOME) WRITE(IW,9992)
C
      CUTOFF=TOL
      NOCA  = NOA-NACORE
      NVIRA = NO-NOA
      NORB  = NO-NACORE
      NOC   = NOA
      NVIR  = NVIRA
C
      NBF2  = (NBF*NBF+NBF)/2
      NBF3  = NBF*NBF
      NORB2 = (NORB*NORB+NORB)/2
      NROT  = NVIR*NOC
C
C        ---- OPTION TO DEFEAT SYMMETRY USAGE ----
C
      ABEL=ABELPT()
      IF(NOSYMT.EQ.1) CALL SYMOFF
C
C     ----- OBTAIN ORBITAL SYMMETRY LABELS -----
C
      CALL VALFM(LOADFM)
      IVEC  = 1     + LOADFM
      IIREP = IVEC  + NBF3
      LAST  = IIREP + NBF
      NEEDA = LAST  - LOADFM
      CALL GETFM(NEEDA)
C
      CALL VALFM(LOADFM)
      IMOLAB = 1      + LOADFM
      IMODEG = IMOLAB + NBF
      IS     = IMODEG + NBF
      IQ     = IS     + NBF2
      IT     = IQ     + NBF3
      LAST   = IT     + NBF
      NEEDB  = LAST   - LOADFM
      CALL GETFM(NEEDB)
      IF(EXETYP.EQ.CHECK) GO TO 90
C
      CALL DAREAD(IDAF,IODA,X(IVEC),NBF3,15,0)
      CALL DAREAD(IDAF,IODA,X(IS  ),NBF2,12,0)
      CALL DAREAD(IDAF,IODA,X(IQ  ),NBF3,45,0)
      CALL TRFSYM(X(IMOLAB),X(IIREP),X(IMODEG),X(IQ),X(IS),X(IVEC),
     *            X(IT),IA,NQMT,NBF,NO,NBF)
C
   90 CONTINUE
      CALL RETFM(NEEDB)
C
      NTSAVE = NT
      IF(DIRSCF  .AND.  .NOT.ABEL) NT=1
C
C     ----- COMPLETE MP2 LAGRANGIAN -----
C
C     -IVEC-   NBF*NBF                       - ORBITAL COEFF     -
C     -IIREP-  NBF                           - MOSIRP            -
C     -INROT-  NROT                          - XLAI1,            -
C     -IRHO-   MAX(NOC*NOC,NVIR*NVIR,NBF2)   - PIJ,PAB,PMN,XLAI2 -
C     -IXLMN-  NBF3                          - X,XLMN            -
C     -IXX-    MAX(NBF3,NINTMX)              - XX,Y              -
C     -IIX-    NINTMX                        - IX                -
C
      CALL VALFM(LOADFM)
      ILAI  = 1     + LOADFM
      IENG  = ILAI  + MAX(NINTMX,NROT)
      LAST  = IENG  + NBF
      NEEDC = LAST  - LOADFM
      CALL GETFM(NEEDC)
C
      CALL VALFM(LOADFM)
      IRHO  = 1     + LOADFM
      IXLMN = IRHO  + MAX(NOC*NOC,NVIR*NVIR,NBF2)
      IXX   = IXLMN + NBF3
      IIX   = IXX   + MAX(NBF3,NINTMX)
      LAST  = IIX   + NINTMX
      NEEDD = LAST  - LOADFM
      CALL GETFM(NEEDD)
      CALL TSECND(TIM0)
      IF(EXETYP.EQ.CHECK) GO TO 250
C
      CALL DAREAD(IDAF,IODA,X(IENG),NBF,17,0)
      CALL DAREAD(IDAF,IODA,X(IVEC),NBF3,15,0)
C
      IDEN = IRHO
      IFCK = IXLMN
      IBUF = IXX
      CALL MP2LAI(X(IRHO),X(IRHO),X(IRHO),X(IXLMN),X(IXX),X(IXLMN),
     *            X(IXX),X(IIX),X(ILAI),X(IRHO),NOC,NVIR,NBF,
     *            KREC1,KREC2,KREC3,KREC4,X(IVEC))
C
      IF(DBUG) THEN
         IF(SOME) WRITE(IW,9999)
         CALL PRSQ(X(ILAI),NOC,NVIR,NVIR)
      END IF
C
  250 CONTINUE
      CALL TSECND(TIM1)
      TIME = TIM1 - TIM0
      IF(SOME) WRITE(IW,9997) TIME
      CALL FLSHBF(IW)
C
      CALL RETFM(NEEDD)
C
C     ----- MP2 CPHF EQUATIONS -----
C
      IF(MOCPHF.EQ.1) THEN
C
         CALL VALFM(LOADFM)
         IXX   = 1    + LOADFM
         IIX   = IXX  + MAX(NROT,NORB2)
         IT    = IIX  + MAX(NROT,NORB2)
         IB    = IT   + 5*NROT
         IWXY  = IB   + NROT
         LAST  = IWXY + NROT*NROT
         NEED1 = IWXY - LOADFM
         NEEDE = LAST - LOADFM
C
         CALL GOTFM(NGOTMX)
         IF(NWORD.NE.0) NGOTMX = MIN(NWORD - NEEDA - NEEDC,NGOTMX)
C
C     COMPILATION OF THE NEXT LINE FORCES SMALL PROBLEMS INTO A
C     TWO SLICE, OUT-OF-CORE ORBITAL HESSIAN STORAGE ALGORITHM.
C
C-DBG-    NGOTMX = NEED1 + (NROT*NROT+10)/2
C
         CALL DAREAD(IDAF,IODA,X(IENG),NBF,17,0)
C
         IF(NEEDE.GT.NGOTMX) THEN
            MEM = NGOTMX - NEED1
            IF(MEM.LT.NROT) THEN
               IF(MASWRK) WRITE(IW,9998) NROT-MEM
               CALL ABRT
            END IF
            NEEDE  = NEED1 + MEM
            CALL GETFM(NEEDE)
            IF(EXETYP.EQ.CHECK) GO TO 350
C
            J00=IWXY
            J01=IIX
            J02=IB-NROT
C
            CALL MPCPHO(X(IWXY),X(IB),X(IXX),X(IIX),X(IENG),
     *                  X(IIREP),X(ILAI),NORB,NORB2,NOC,NVIR,
     *                  X(IT),X(IB),NROT,NFTIJ,NFTIA,CUTOFF,MEM)
         ELSE
            CALL GETFM(NEEDE)
            IF(EXETYP.EQ.CHECK) GO TO 350
            CALL MPCPHF(X(IWXY),X(IWXY),X(IXX),X(IIX),X(IENG),
     *                  X(IIREP),X(ILAI),NORB,NORB2,NOC,NVIR,
     *                  X(IT),X(IB),NROT,NFTIJ,NFTIA)
         END IF
C
      ELSE
C
         WRITE(IW,9980)
C
         MAXIO  = 2*MAXDII
         ITERMX = MAXITC
         ITERMX2= (ITERMX*ITERMX+ITERMX)/2
C
         CALL VALFM(LOADFM)
         IRHO  = 1     + LOADFM
         IXX   = IRHO  + NBF2
         IY    = IXX   + NBF3
         IXLMN = IY    + NBF3
         LBUF  = IXLMN + NBF3
         IBUF  = LBUF  + NINTMX
         IPAI  = IBUF  + NINTMX
         IPBJ  = IPAI  + NROT
         IERR  = IPBJ  + NROT
         IWRK  = IERR  + NROT
         IAM   = IWRK  + NROT
         IX    = IAM   + MAXDII*MAXDII
         IPVT  = IX    + ITERMX
         IB    = IPVT  + ITERMX
         IODIIS= IB    + ITERMX2
         LAST  = IODIIS+ MAXIO
         NEEDE = LAST  - LOADFM
         CALL GETFM(NEEDE)
         IF(EXETYP.EQ.CHECK) GO TO 350
C
         CALL DAREAD(IDAF,IODA,X(IENG),NBF,17,0)
C
         CALL AOCPHF(X(IRHO),X(IXX),X(IY),X(IXLMN),X(LBUF),X(IBUF),
     *               X(ILAI),X(IPAI),X(IPBJ),NOC,NVIR,NBF,X(IVEC),
     *               X(IERR),X(IWRK),X(IAM),X(IX),X(IPVT),X(IB),
     *               X(IODIIS),X(IENG),NROT,ITERMX)
C
      CALL DCOPY(NROT,X(IPAI),1,X(ILAI),1)
C
      END IF
C
  350 CONTINUE
      CALL RETFM(NEEDE)
C
      NT = NTSAVE
      IF(NOSYMT.EQ.1) CALL SYMON
C
      CALL TSECND(TIM2)
      TIME = TIM2 - TIM1
      IF(SOME) WRITE(IW,9996) TIME
      CALL FLSHBF(IW)
C
C     ----- COMPLETE W2 -----
C
C     ILAI  =  NROT         PAI, OBTAINED FROM CPHF EQUATION
C     IENG  =  NBF          ORBITAL ENERGY
C     IPOCC =  NOC*NOC      OCCUPIED-OCCUPIED DENSITY (POCC)
C     IPVIR =  NVIR*NVIR    VIRTUAL-VIRTUAL DENSITY (PVIR)
C     IWOCC =  NOC*NOC      ENERGY WEIGHT DENSITY (WOCC)
C     IWAI  =  NOC*NVIR     ENERGY WEIGHT DENSITY (WAI)
C     IXX   =  NBF2         XX
C     IX    =  NBF3         X
C     IWRK1 =  NBF3         SCRATCH
C     IWRK2 =  NBF3         SCRATCH
C
      CALL VALFM(LOADFM)
      IPOCC   = 1   + LOADFM
      IPVIR   = IPOCC + NOC*NOC
      IWOCC   = IPVIR + NVIR*NVIR
      IWAI    = IWOCC + NOC*NOC
      IXX     = IWAI  + NVIR*NOC
      IX      = IXX   + NBF2
      IWRK1   = IX    + NBF3
      IWRK2   = IWRK1 + NBF3
      LAST    = IWRK2 + NBF3
      IF(MOCPHF.EQ.0) THEN
         NBUF    = LAST
         IBUF    = NBUF + NINTMX
         LAST    = IBUF + NINTMX
      END IF
      NEEDF   = LAST  - LOADFM
      CALL GETFM(NEEDF)
      IF(EXETYP.EQ.CHECK) GO TO 550
C
      CALL DAREAD(IDAF,IODA,X(IPOCC),NOC *NOC ,KREC1,0)
      CALL DAREAD(IDAF,IODA,X(IPVIR),NVIR*NVIR,KREC2,0)
      CALL DAREAD(IDAF,IODA,X(IWAI ),NVIR*NOC ,KREC3,0)
      CALL DAREAD(IDAF,IODA,X(IWOCC),NOC *NOC ,KREC5,0)
C
      IF(MOCPHF.EQ.1) THEN
         CALL MPW2IJ(X(IPOCC),X(IPVIR),X(ILAI),X(IWOCC),X(IWAI),
     *               X(IENG),X(IXX),X(IX),NOC,NVIR,NORB,
     *               NORB2,NFTIJ,NFTIA)
      END IF
C
C     ----- TRANSFORM PIJ, PAB, AND PAI TO P2(MUNU) -----
C
      CALL MPGTRP(X(IPOCC),X(IPVIR),X(ILAI),X(IVEC),NBF,NOC,NVIR,
     *            X(IWRK1),X(IWRK2),X(IXX))
      CALL DAWRIT(IDAF,IODA,X(IXX),NBF2,KREC7,0)
C
      IF(MOCPHF.EQ.0) THEN
         CALL W2IJAO(X(IPOCC),X(ILAI),X(IWOCC),X(IWAI),X(IENG),
     *               X(IXX),X(IX),X(IWRK1),X(IWRK2),X(NBUF),X(IBUF),
     *               X(IVEC),NOC,NVIR,NBF,NBF2)
      END IF
C
      CALL DAREAD(IDAF,IODA,X(IX),NVIR*NVIR,KREC6,0)
C
C     ----- TRANSFORM WIJ, WAB, AND WAI TO P2(MUNU) -----
C
      CALL MPGTRP(X(IWOCC),X(IX),X(IWAI),X(IVEC),NBF,NOC,NVIR,
     *            X(IWRK1),X(IWRK2),X(IXX))
      CALL MPWMP2(X(IVEC),X(IENG),X(IWRK1),X(IXX),NOC,NBF,
     *            KREC7,KREC8,KREC9)
C
      CALL TSECND(TIM3)
      TIME = TIM3 - TIM2
      IF(SOME) WRITE(IW,9995) TIME
      IF(MASWRK)
     *    WRITE(IW,*) '..... DONE WITH MP2 CPHF AND DENSITY .....'
      CALL TIMIT(1)
      CALL FLSHBF(IW)
C
  550 CONTINUE
      CALL RETFM(NEEDF)
C
      IF(NPRINT.NE.-5) CALL MP2NOS
C
      ETOT = EMP2
      CALL DAWRIT(IDAF,IODA,ENUC,MXRT+15,2,0)
C
C     COSMO CHANGES - KIM, PUT IN BY VOLKER NOV.1998
C
C     COMPUTE PROPERTIES INCLUDING THE COSMO SOLVATION CORRECTION
C     FOR THE SECOND ORDER WAVE FUNCTION AND ITERATE BOTH Q'S
C     AND MP2 ENERGY TO SELF-CONSISTENCY
C
      IF(ISEPS) THEN
          WRITE(IW,*) "MP2ITER IN MP2GRD = ",MP2ITER
          WRITE(IW,*)
          WRITE(IW,*)"AMOS'S DISTRIBUTED MULTIPOLE ANALYSIS:"
          WRITE(IW,*)"      -- CALLING FROM MP2GRD --"
          WRITE(IW,*)
C
          CALL NUMPROP
C
          WRITE(IW,*)"NEW SURFACE CHARGES/SCREENING E FOR MP2",MP2TRIP
          IF(MP2TRIP.EQ.1) THEN
             NDEN2=(NDEN*NDEN+NDEN)/2
             LENAB2 = NPS2 + NDEN*NPS + NDEN2
             CALL VALFM(LOADFM)
             LABCMT = LOADFM + 1
             LAST   = LABCMT + LENAB2
             NEEDCO = LAST - LOADFM - 1
             CALL GETFM(NEEDCO)
C
C              RECOVER PREVIOUS A,B,C MATRICES
C
             NFT26 = 26
             LENB = NDEN*NPS
             CALL SQREAD(NFT26,X(LABCMT)              ,NPS2)
             CALL SQREAD(NFT26,X(LABCMT+NPS2)         ,LENB)
             CALL SQREAD(NFT26,X(LABCMT+NPS2+NDEN*NPS),NDEN2)
             CALL SEQREW(NFT26)
C
C              AND COMPUTE THE SOLVATION ENERGY WITH THEM
C
             CALL DECORR(X(LABCMT),LENAB2)
             CALL RETFM(NEEDCO)
C
             DSECORR=SECORR-SAVESE
C            EMP2SOL1=(EMP2-ETOTSMP)*627.517D+00
C            EMP2SOL=(EMP2-ETOTSMP-SECORR)*627.517D+00
C            EMP2SOL2=(EMP2-ETOTSMP-DSECORR)*627.517D+00
      CALL FLSHBF(6)
C
C   NOW WE HAVE TO DO THE OCE CORRECTION ANALOGOUS TO THE SCF PART
C   QVCOSMO SHOULD BE THERE FROM DECORR, WE NEED QVGAMESS FROM THE
C   MP2 DENSITY - THIS IS FOR THE CONVERGED MP2 RESULTS, MP2ITER=0
C
             EOC1=0.0D+00
             IF(MP2ITER.EQ.0) THEN
                DO 5 II=1,NPS
                   QVPOT(II)=0.0D+00
  5             CONTINUE
      CALL FLSHBF(6)
C     CALL VALFM(LOADFM)
C     LAST  = 1   + LOADFM
C     LAST  = LAST + NQS
C     NEEDG2= LAST - LOADFM
C     CALL GETFM(NEEDG2)
                CALL COSPOT(QVPOT)
C
C   COSPOT GAVE US V'; NOW FORM QV'
C
                QVCORR=0.0D+00
                DO 6 I=1,NQS
                   QVCORR = QVCORR + QVPOT(I)*COSZAN(I)
  6             CONTINUE
                WRITE(IW,*)"QV (GAMESS,MP2): QVCORR = ",QVCORR
                WRITE(IW,*)"QV (COSMO,MP2) : QVCOSMO= ",QVCOSMO
                EOC1=-QVCORR+QVCOSMO
                EOC2=EOC1*627.517D+00
                CHECK2=QVCORR*QVCORR-QVCOSMO*QVCOSMO
                WRITE(IW,*)"-QV(GAMESS,MP2)+QV(COSMO,MP2):",EOC1
                WRITE(IW,*)"POTENTIAL NORM DIFFERENCE :",CHECK2
                WRITE(IW,*)"QVCOSMO/QVGAMESS:",QVCOSMO/QVCORR
                WRITE(IW,*)
                WRITE(IW,9400)
                WRITE(IW,*)"  OUTLYING CHARGE EFFECT (A.U.)      =",EOC1
                WRITE(IW,*)"  OUTLYING CHARGE EFFECT (KCAL/MOL)  =",EOC2
                WRITE(IW,9400)
             ENDIF
      CALL FLSHBF(6)
C
C   EOC1 IS NOW THE MP2 OUTLYING CHARGE CORRECTION, WHICH MUST BE
C   ADDED TO EMP2 AFTER PRINTOUT
C
             EMP2C=EMP2+EOC1
             WRITE(IW,9500)ESCF,SAVESE,SECORR,DSECORR,EMP2,EMP2C
             EMP2=EMP2C
         ENDIF
C
         IF(MP2ITER.EQ.1) THEN
             WRITE(IW,*)" 1 ELECTRON CONTRIBUTION FROM MP2 CHARGES:"
             CALL ONEEI
         ENDIF
C
         WRITE(IW,*)
         WRITE(IW,*) "**********COSMO INFORMATION**MP2GRD**************"
         WRITE(IW,*)
         WRITE(IW,*)"  NUMBER OF SURFACE SEGMENTS IS ",NPS
         WRITE(IW,*)
         AREA=0.0D+00
         DO 889 I=1,NPS
           AREA=AREA+AR(I)
           WRITE(IW,'(2I4,6F10.5)') I,IATSP(I),CORZAN(1,I),CORZAN(2,I),
     *CORZAN(3,I),COSZAN(I),AR(I),COSZAN(I)/AR(I)
 889       CONTINUE
           WRITE(IW,*)
           WRITE(IW,*)" TOTAL SURFACE AREA OF CAVITY (A**2)     = ",AREA
           WRITE(IW,*)
           WRITE(IW,*)"***********************************************"
           WRITE(IW,*)
C
C     MP2ITER=1 MEANS THAT THE MP2 ENERGY IS NOT CONVERGED, SO WE GO
C     BACK HERE AFTER THE MP2 DENSITY IS UPDATED
C
         IF(MP2ITER.EQ.1) THEN
            WRITE(IW,*) "RETURNING FROM MP2GRD, MP2ITER=",MP2ITER
            CALL RETFM(NEEDC)
            CALL RETFM(NEEDA)
            RETURN
         END IF
      END IF   ! END OF COSMO
C
C        A RUN TO GENERATE ONLY THE DENSITY CAN NOW EXIT.
C
      IF(MPPROP.EQ.1) THEN
         CALL RETFM(NEEDC)
         CALL RETFM(NEEDA)
         RETURN
      END IF
C
C     ----- NON-SEPARABLE CONTRIBUTIONS TO THE 2ND ORDER DENSITY -----
C     IAIBJ  = NVIR*NOC*NVIR       AIBJ
C     IWRK1  = NBF3                X
C     IWRK2  = NBF3                Y,XX
C     IXIN   = NBF3                XIN
C
      CALL VALFM(LOADFM)
      IAIBJ = 1   + LOADFM
      IWRK1 = IAIBJ + NVIR*NOCA*NVIR
      IWRK2 = IWRK1 + NBF3
      IXIN  = IWRK2 + NBF3
      IENG  = IXIN  + NBF2*NBF
      LAST  = IENG  + NBF
      NEEDG = LAST - LOADFM
C
      CALL GOTFM(NGOTMX)
      IF(NEEDG.GT.NGOTMX) THEN
         IF(MASWRK) WRITE(IW,9993) NEEDG-NGOTMX
         CALL ABRT
      END IF
      CALL GETFM(NEEDG)
      IF(EXETYP.EQ.CHECK) GO TO 750
C
      CALL DAREAD(IDAF,IODA,X(IENG),NBF,17,0)
      CALL MAMNLS(X(IAIBJ),X(IENG),X(IVEC),X(IWRK1),X(IWRK2),X(IWRK2),
     *            X(IXIN),NBF,NBF2,NOC,NVIR,NFTIA,NFT16,
     *            CUTOFF,SOME,NACORE,NOCA)
C
  750 CONTINUE
      CALL RETFM(NEEDG)
      CALL RETFM(NEEDC)
      CALL RETFM(NEEDA)
C
      IF(MASWRK) WRITE(IW,*) '..... DONE WITH MP2 2DM .....'
      CALL TIMIT(1)
      RETURN
C
 9999 FORMAT(/' --- LAGRANGIAN ---'/)
 9998 FORMAT(/' OUT OF MEMORY FOR CPHF EQUATIONS'/
     *        ' INCREASE MEMORY SIZE BY AT LEAST ',I10,' WORDS'/)
 9997 FORMAT(/' -MP2- LAGRANGIAN     TOOK = ',F10.3,' SECONDS')
 9996 FORMAT( ' -MP2- CPHF EQUATIONS TOOK = ',F10.3,' SECONDS')
 9995 FORMAT( ' -MP2- W2(OCC,OCC)    TOOK = ',F10.3,' SECONDS')
 9993 FORMAT(/' OUT OF MEMORY FOR 2DM TRANSFORMATION'/
     *        ' INCREASE MEMORY SIZE BY AT LEAST ',I10,' WORDS'/)
 9992 FORMAT(/10X,21(1H-)/10X,'CPHF FOR MP2 GRADIENT',
     *       /10X,21(1H-))
 9980 FORMAT(/' ..... SOLVING CPHF IN THE AO BASIS .....')
 9400 FORMAT(2X,'***************************************************',
     *'**************')
 9500 FORMAT(/1X,'RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE'/
     *       8X,'  ECOSMO(0)                      =',1X,F20.10/
     *       8X,'  RHF SECORR                     =',1X,F20.10/
     *       8X,'  MP2 SECORR                     =',1X,F20.10/
     *       8X,'DELTA SECORR                     =',1X,F20.10/
     *       8X,'E(COSMO-MP2)                     =',1X,F20.10/
     *       8X,'E(COSMO-MP2,OCE-CORRECTED)       =',1X,F20.10/)
      END
C
C*MODULE MP2GRD  *DECK MP2LAI
      SUBROUTINE MP2LAI(PIJ,PAB,PMN,X,Y,XLMN,BUF,IBUF,XLAI1,XLAI2,
     *                  NOC,NVIR,NBF,KREC1,KREC2,KREC3,KREC4,C)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E
C
      DIMENSION PIJ(NOC,NOC), PAB(NVIR,NVIR), PMN(*), XLAI2(NVIR,NOC)
      DIMENSION X(NBF,NBF),   XLMN(*)
      DIMENSION Y(NBF,NBF),   BUF(NINTMX), IBUF(NINTMX)
      DIMENSION XLAI1(NVIR,NOC), C(NBF,*)
C
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- COMPLETE CONSTRUCTION OF THE MP2 LAGRANGIAN -----
C     ACCORDING TO EQUATION 11 OF THE FRISCH PAPER.
C
C        BEGIN TRANSFORMING -PIJ- AND -PAB- TO AO BASIS AS -PMN-
C
      CALL DAREAD(IDAF,IODA,PIJ,NOC*NOC,KREC1,0)
      CALL MRARTR(PIJ,NOC,NOC,NOC,C,NBF,NBF,X,NBF)
      CALL MRARBR(C,NBF,NBF,NOC,X,NBF,NBF,Y,NBF)
C
      CALL DAREAD(IDAF,IODA,PAB,NVIR*NVIR,KREC2,0)
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,
     *           ONE,Y,NBF)
C
      MUNU = 0
      DO 110 MU=1,NBF
         DO 100 NU=1,MU
            MUNU=MUNU+1
            PMN(MUNU) = Y(MU,NU) + Y(NU,MU)
  100    CONTINUE
  110 CONTINUE
C
C     THIRD AND FOURTH TERM LAGRANGIAN CONTRIBUTIONS, IN AO BASIS
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
C     ---- TRANSFORM THESE TWO TERMS TO MO BASIS,
C          AND ADD ALL CONTRIBUTIONS TOGETHER ----
C
      CALL DAREAD(IDAF,IODA,XLAI1,NVIR*NOC,KREC4,0)
C
      CALL MRARBR(Y,NBF,NBF,NBF,C,NBF,NOC,X,NBF)
      CALL DGEMM('T','N',NVIR,NOC,NBF,ONE,C(1,NOC+1),NBF,X,NBF,
     *           ONE,XLAI1,NVIR)
C
      CALL DAREAD(IDAF,IODA,XLAI2,NVIR*NOC,KREC3,0)
C
      DO 30 MI=1,NOC
      DO 30 MA=1,NVIR
   30    XLAI1(MA,MI) = -( XLAI1(MA,MI)+XLAI2(MA,MI) )
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK P2W2N
      SUBROUTINE P2W2N(XPQRJ,T2,TOL,E,POCC,WOCC,PVIR,WVIR,XLAI1,XLAI2,
     *                 NORB,NOC,NVIR,NOCC,NOC1,NBF2,NBF,NCO,NOCA)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
      PARAMETER (ZERO=0.0D+00, PT1=0.1D+00, TWO=2.0D+00)
C
      COMMON /IJPAIR/ IA(MXAO)
C
      DIMENSION XPQRJ(NOCC,NBF2,NBF),T2(NVIR,NOC,NOCC)
      DIMENSION E(NORB),POCC(NOC,NOC),WOCC(NOC,NOC),
     *          PVIR(NVIR,NVIR),WVIR(NVIR,NVIR),
     *          XLAI1(NVIR,NOC),XLAI2(NVIR,NOC)
C
      TOL1 = TOL*PT1
C
      DO 150 MB=1,NVIR
C
C     ----- FIRST THE AMPLITUDES ------
C
      DO 104 MA=1,NVIR
         MAJ = IA(MA+NOC)
         MBJ = IA(MB+NOC)
         DO 102 MJ=1,NOC
            MAJ = MAJ+1
            MBJ = MBJ+1
            DO 100 MK=1,NOCC
               DJKAB = E(MJ) + E(MK+NCO+NOC1-1) -
     *                 E(MA+NOC) - E(MB+NOC)
               TERM = TWO*XPQRJ(MK,MAJ,MB+NOC) -
     *                XPQRJ(MK,MBJ,MA+NOC)
               TERM = TERM/DJKAB
               IF(ABS(TERM).LT.TOL1) TERM = ZERO
               T2(MA,MJ,MK) = TERM
  100       CONTINUE
  102    CONTINUE
  104 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO POCC, WOCC, AND LAI(II)
C
      DO 120 MJ=NCO+1,NOC
         DO 118 MK=1,NOCC
            DO 114 MA=1,NVIR
C
               TERM = T2(MA,MJ,MK)
               IF(TERM.EQ.ZERO) GO TO 114
C
               MAI = IA(MA+NOC)
               DO 110 MI=1,NOC
                  MAI = MAI+1
                  MIJ = IA(MAX(MI,MJ)) + MIN(MI,MJ)
                  DIKAB = E(MI) + E(MK+NCO+NOC1-1) - E(MA+NOC) -
     *                    E(MB+NOC)
                  XLAI2(MA,MI) = XLAI2(MA,MI) -
     *                           TERM*XPQRJ(MK,MIJ,MB+NOC)
                  WOCC(MI,MJ) = WOCC(MI,MJ) -
     *                          TERM*XPQRJ(MK,MAI,MB+NOC)
                  IF(MI.LE.NCO) THEN
                     DIJ = E(MJ) - E(MI)
                     POCC(MI,MJ) = POCC(MI,MJ) +
     *                             TERM*XPQRJ(MK,MAI,MB+NOC)/DIJ
                  ELSE
                     POCC(MI,MJ) = POCC(MI,MJ) -
     *                             TERM*XPQRJ(MK,MAI,MB+NOC)/DIKAB
                  END IF
  110          CONTINUE
  114       CONTINUE
  118    CONTINUE
  120 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO PVIR, WVIR, AND LAI(I)
C
      DO 140 MI=1,NOCA
         DO 138 MK=1,NOCC
            DO 134 MC=1,NVIR
C
               TERM = T2(MC,MI+NCO,MK)
               IF(TERM.EQ.ZERO) GO TO 134
C
               DO 130 MA=1,NVIR
                  MAI = IA(MA+NOC) + MI+NCO
                  MAC = IA(MAX(MA,MC)+NOC) + MIN(MA,MC) + NOC
                  DIKAB = E(MI+NCO) + E(MK+NCO+NOC1-1) -
     *                    E(MA+NOC) - E(MB+NOC)
                  PVIR(MA,MC) = PVIR(MA,MC) +
     *                          TERM*XPQRJ(MK,MAI,MB+NOC)/DIKAB
                  WVIR(MA,MC) = WVIR(MA,MC) -
     *                          TERM*XPQRJ(MK,MAI,MB+NOC)
                  XLAI1(MA,MI+NCO)= XLAI1(MA,MI+NCO) +
     *                          TERM*XPQRJ(MK,MAC,MB+NOC)
  130          CONTINUE
  134       CONTINUE
  138    CONTINUE
  140 CONTINUE
C
  150 CONTINUE
C
      IF(NCO.GT.0) THEN
         DO 200 MJ=1,NCO
            DO 210 MI=1,NOC
                POCC(MI,MJ) = POCC(MJ,MI)
  210       CONTINUE
  200    CONTINUE
      END IF
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK SPARSQ
      SUBROUTINE SPARSQ(A,LDA,B,X,N,G,H,XI,XJ,SPARAX,SPARXA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION A(LDA,*),X(N),B(N),G(N),H(N),XI(N),XJ(N)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      EXTERNAL SPARAX,SPARXA
C
      DATA ZERO /0.0D+00/, EPS/1.0D-07/
C
C     ----- SOLVE A * X = B , USING AN ITERATIVE PROCEDURE       -----
C
C     ----- NUMERICAL RECIPES (P.70), CAMBRIDGE UNIVERSITY PRESS -----
C          W.H.PRESS, B.P.FLANNERY, S.A.TEUKOLSKY, W.T.VETTERLING
C
      EPS2=N*EPS*EPS
      IRST=0
   10 IRST=IRST+1
      CALL SPARAX(A,LDA,X,XI,N)
      RP=ZERO
      BSQ=ZERO
      DO 11 J=1,N
         BSQ=BSQ+B(J)*B(J)
         XI(J)=XI(J)-B(J)
         RP=RP+XI(J)*XI(J)
   11 CONTINUE
      CALL SPARXA(A,LDA,XI,G,N)
      DO 12 J=1,N
      G(J)=-G(J)
      H(J)= G(J)
   12 CONTINUE
      DO 19 ITER=1,10*N
         CALL SPARAX(A,LDA,H,XI,N)
         ANUM=ZERO
         ADEN=ZERO
         DO 13 J=1,N
            ANUM=ANUM+G(J)*H(J)
            ADEN=ADEN+XI(J)*XI(J)
   13    CONTINUE
         IF(ADEN.EQ.ZERO) THEN
           IF(MASWRK) WRITE(IW,9010)
           CALL ABRT
         END IF
         ANUM=ANUM/ADEN
         DO 14 J=1,N
            XI(J)=X(J)
            X(J)=X(J)+ANUM*H(J)
   14    CONTINUE
         CALL SPARAX(A,LDA,X,XJ,N)
         RSQ=ZERO
         DO 15 J=1,N
            XJ(J)=XJ(J)-B(J)
            RSQ=RSQ+XJ(J)*XJ(J)
   15    CONTINUE
         IF(RSQ.EQ.RP.OR.RSQ.LE.BSQ*EPS2) RETURN
         IF(RSQ.GT.RP) THEN
            DO 16 J=1,N
               X(J)=XI(J)
   16       CONTINUE
            IF(IRST.GE.3) RETURN
            GO TO 10
         END IF
         RP=RSQ
         CALL SPARXA(A,LDA,XJ,XI,N)
         GG=ZERO
         DGG=ZERO
         DO 17 J=1,N
            GG=GG+G(J)*G(J)
            DGG=DGG+(XI(J)+G(J))*XI(J)
   17    CONTINUE
         IF(GG.EQ.ZERO) RETURN
         GAM=DGG/GG
         DO 18 J=1,N
            G(J)=-XI(J)
            H(J)=G(J)+GAM*H(J)
   18    CONTINUE
   19 CONTINUE
C
      IF(MASWRK) WRITE(IW,9020)
      CALL ABRT
      RETURN
C
 9010 FORMAT(1X,'VERY SINGULAR MATRIX IN SPARSQ')
 9020 FORMAT(1X,'TOO MANY ITERATIONS IN SPARSQ')
      END
C
C*MODULE MP2GRD  *DECK SPARXA1
      SUBROUTINE SPARXA1(A,LDA,X,XA,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      DIMENSION A(LDA,*),X(N),XA(N)
      CALL DGEMV('T',N,N,ONE,A,LDA,X,1,ZERO,XA,1)
      RETURN
      END
C
C*MODULE MP2GRD  *DECK SPARAX1
      SUBROUTINE SPARAX1(A,LDA,X,AX,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      DIMENSION A(LDA,*),X(N),AX(N)
      CALL DGEMV('N',N,N,ONE,A,LDA,X,1,ZERO,AX,1)
      RETURN
      END
C
C*MODULE MP2GRD  *DECK TRFOUT
      SUBROUTINE TRFOUT(X,NX,IX,NIX,N,I,J)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),IX(*)
C
      PARAMETER (NOMAX=3)
      DIMENSION IO(NOMAX),JO(NOMAX),KO(NOMAX),LO(NOMAX),
     *          IJKLO(NOMAX),OIJKL(NOMAX)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PCKLAB/ LABSIZ
C
C     ----- PRINT ORDERED -MO- INTEGRALS AS WRITTEN OUT -----
C
      WRITE(IW,9999) I,J
      IF(NX.EQ.0) GO TO 60
C
*I64  MASK32 = 2**32 - 1
C
      MO=NIX
      NO=NIX*2
      IF(MOD(NX,2).EQ.1) THEN
         NO=NO-1
         IX(NO)=IX(MO)
         NO=NO-1
         MO=MO-1
      END IF
C
      IF (LABSIZ .EQ. 2) THEN
*I32     MO = MO/2
*I64     MO = MO
      ELSE IF (LABSIZ .EQ. 1) THEN
*I32     MO = MO
*I64     IF (MOD(NO/2,2) .EQ. 0) THEN
*I64       MO = MO*2
*I64     ELSE
*I64       MO = (MO*2) - 1
*I64     END IF
      END IF
C
      DO 10 M = MO, 1, -1
C
                     NPACK = M
                     IF (LABSIZ .EQ. 2) THEN
*I32                   IPACK = IX( 2*NPACK - 1 )
*I32                   JPACK = IX( 2*NPACK     )
*I64                   LABEL = IX(NPACK)
*I64                   IPACK = ISHFT( LABEL, -32 )
*I64                   JPACK = IAND( LABEL, MASK32 )
                     ELSE IF (LABSIZ .EQ. 1) THEN
*I32                   LABEL = IX(NPACK)
*I32                   IPACK = ISHFT( LABEL, -16 )
*I32                   JPACK = IAND( LABEL, 65535 )
*I64                   IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                     LABEL = IX( NPACK/2 )
*I64                     IPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     JPACK = IAND( LABEL, 65535 )
*I64                   ELSE
*I64                     LABEL = IX( (NPACK/2)+1 )
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                   END IF
                     END IF
                     IWORD2 = IPACK
                     IWORD1 = JPACK
C
         IX( 2*M-1)=IWORD2
         IX( 2*M  )=IWORD1
   10 CONTINUE
C
C     ----- PRINT OUT -----
C
      MO=1
      NO=1
      KL=0
      DO 50 K=1,N
      DO 40 L=1,K
      KL=KL+1
      IF(KL.NE.IX(MO)) GO TO 40
      IO(NO)=I
      JO(NO)=J
      KO(NO)=K
      LO(NO)=L
      IJKLO(NO)=  KL
      OIJKL(NO)=X(MO)
      MO=MO+1
      NO=NO+1
      IF(NO.LE.NOMAX) GO TO 40
      NO=NO-1
      WRITE(IW,9998) (IO(M),JO(M),KO(M),LO(M),IJKLO(M),
     *                        OIJKL(M),M=1,NO)
      NO=1
   40 CONTINUE
   50 CONTINUE
      IF(NO.EQ.1) GO TO 60
      NO=NO-1
      WRITE(IW,9998) (IO(M),JO(M),KO(M),LO(M),IJKLO(M),
     *                        OIJKL(M),M=1,NO)
C
   60 CONTINUE
      RETURN
 9999 FORMAT(' --- I,J = ',2I5)
 9998 FORMAT(3(4I3,6X,I5,F20.12))
      END
C*MODULE MP2GRD  *DECK TRFPRT
      SUBROUTINE TRFPRT(X,IX,N,I,J)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),IX(*)
      DIMENSION IO(3),JO(3),KO(3),LO(3),IJKLO(3),OIJKL(3)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      DATA NOMAX /3/
C
C     ----- PRINT ORDERED -AO- INTEGRALS AS READ IN -----
C
      WRITE(IW,9999) I,J
C
      M=1
      NO=1
      KL=0
      DO 15 K=1,N
         DO 10 L=1,K
            KL=KL+1
            IF(KL.NE.IX(M)) GO TO 10
            M=M+1
            IO(NO)=I
            JO(NO)=J
            KO(NO)=K
            LO(NO)=L
            IJKLO(NO)=KL
            OIJKL(NO)=X(KL)
            NO=NO+1
            IF(NO.LE.NOMAX) GO TO 10
            NO=NO-1
            WRITE(IW,9998) (IO(ID),JO(ID),KO(ID),LO(ID),
     *                     IJKLO(ID),OIJKL(ID),ID=1,NO)
            NO=1
   10    CONTINUE
   15 CONTINUE
C
      IF(NO.EQ.1) GO TO 20
      NO=NO-1
      WRITE(IW,9998) (IO(ID),JO(ID),KO(ID),LO(ID),
     *                IJKLO(ID),OIJKL(ID),ID=1,NO)
   20 CONTINUE
      RETURN
C
 9999 FORMAT(' --- I,J = ',2I5)
 9998 FORMAT(3(4I3,6X,I5,E20.12))
      END
C
C*MODULE MP2GRD  *DECK TRFRD
      SUBROUTINE TRFRD(NFT,X,NX,IX,NIX,NN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),IX(*)
C
      COMMON /PCKLAB/ LABSIZ
C
      PARAMETER (ZERO=0.0D+00)
C
      DO 10 I=1,NN
         IX(I)=0
          X(I)=ZERO
   10 CONTINUE
C
*I64  MASK32 = 2**32 - 1
C
      READ(NFT) NX,NIX
      IF(NX.EQ.0) GO TO 60
      CALL TRFSQR(NFT,X,NX,IX,NIX)
C
      MX =NX
      MIX=NIX
C
      IF(MOD(NX,2).EQ.1) THEN
         IX(MX)=IX(MIX)
         MIX= MIX-1
         MX = MX -1
      END IF
C
      IF (LABSIZ .EQ. 2) THEN
*I32     M = MIX/2
*I64     M = MIX
      ELSE IF (LABSIZ .EQ. 1) THEN
*I32     M = MIX
*I64     IF (MOD(MX,4) .EQ. 0) THEN
*I64        M = MIX*2
*I64     ELSE
*I64        M = MIX*2 - 1
*I64     END IF
      END IF
      DO 20 I = M, 1, -1
C
         NPACK = I
         IF (LABSIZ .EQ. 2) THEN
*I32        IPACK = IX( 2*NPACK - 1 )
*I32        JPACK = IX( 2*NPACK     )
*I64        LABEL = IX(NPACK)
*I64        IPACK = ISHFT( LABEL, -32 )
*I64        JPACK = IAND( LABEL, MASK32 )
         ELSE IF (LABSIZ .EQ. 1) THEN
*I32        LABEL = IX(NPACK)
*I32        IPACK = ISHFT( LABEL, -16 )
*I32        JPACK = IAND( LABEL, 65535 )
*I64        IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64           LABEL = IX( NPACK/2 )
*I64           IPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64           JPACK = IAND( LABEL, 65535 )
*I64        ELSE
*I64           LABEL = IX( (NPACK/2)+1 )
*I64           IPACK = ISHFT( LABEL, -48 )
*I64           JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64        END IF
         END IF
C
         IX( 2*I-1 )=IPACK
         IX( 2*I   )=JPACK
   20 CONTINUE
C
C     ----- RECONSTRUCT ROW OF INTEGRALS -----
C
      M=NX
      DO 50 J=1,NN
         I=NN+1-J
         IF(I.NE.IX(M)) X(I)=ZERO
         IF(I.NE.IX(M)) GO TO 50
         X(I)=X(M)
         M=M-1
         IF(M.EQ.0) THEN
            KMAX = I-1
            DO 53 K=1,KMAX
               X(K) = ZERO
   53       CONTINUE
            GO TO 60
         END IF
   50 CONTINUE
C
   60 CONTINUE
      RETURN
      END
C
C*MODULE MP2GRD  *DECK TRFSQR
      SUBROUTINE TRFSQR(NFT,X,NX,IX,NIX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(NX),IX(NIX)
      READ(NFT) X,IX
      RETURN
      END
C
C*MODULE MP2GRD  *DECK TRFSQW
      SUBROUTINE TRFSQW(NFT,X,NX,IX,NIX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NX),IX(NIX)
      WRITE(NFT) X,IX
      RETURN
      END
C
C*MODULE MP2GRD  *DECK TRFWTM
      SUBROUTINE TRFWTM(NFT,X,NX,IX,NIX,NN,I0,TOL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),IX(*)
C
      COMMON/PCKLAB/LABSIZ
C
C     WRITE A ROW OF NON-ZERO INTEGRALS TO DISK FILE -NFT-
C
C     ----- SIFT THE INTEGRALS FOR NON-SMALL ENTRIES -----
C
      NX=0
      DO 10 I=1,NN
         IF( ABS(X(I)).LT.TOL) GO TO 10
         NX    = NX + 1
         X(NX) = X(I)
         IX(NX)= I + I0 - 1
   10 CONTINUE
C
      NIX=NX
      IF(NX.EQ.0) GO TO 40
C
      MX=NX-MOD(NX,2)
      NIX=0
      DO 20 I=1,MX/2
C
                 NPACK = I
                 IPACK = IX( 2*NPACK-1 )
                 JPACK = IX( 2*NPACK   )
                 IF (LABSIZ .EQ. 2) THEN
*I32               IX( 2*NPACK-1 ) = IPACK
*I32               IX( 2*NPACK   ) = JPACK
*I32               NIX = NIX + 2
*I64               LABEL = ISHFT( IPACK, 32 ) + JPACK
*I64               IX(NPACK) = LABEL
*I64               NIX = NIX + 1
                 ELSE IF (LABSIZ .EQ. 1) THEN
*I32               LABEL = ISHFT( IPACK, 16 ) + JPACK
*I32               IX(NPACK) = LABEL
*I32               NIX = NIX + 1
*I64               IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                 LABEL = ISHFT( IPACK, 16 ) + JPACK
*I64                 IX( NPACK/2 ) = IX( NPACK/2 ) + LABEL
*I64               ELSE
*I64                 LABEL = ISHFT( IPACK, 48 ) + ISHFT( JPACK, 32 )
*I64                 IX( (NPACK/2)+1 ) = LABEL
*I64                 NIX = NIX + 1
*I64               END IF
                 END IF
C
   20 CONTINUE
C
      IF(MX.EQ.NX) GO TO 30
      NIX=NIX+1
      IX(NIX)=IX(NX)
C
C     ----- WRITE OUT -----
C
   30 CONTINUE
      WRITE(NFT) NX,NIX
      CALL TRFSQW(NFT,X,NX,IX,NIX)
      RETURN
C
   40 CONTINUE
      WRITE(NFT) NX,NIX
      RETURN
      END
C*MODULE MP2GRD  *DECK W2IJAO
      SUBROUTINE W2IJAO(WOCC3,PAI,WOCC,WAI,E,PMN,XLMN,X,Y,BUF,IBUF,
     *              C,NOC,NVIR,NBF,NBF2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00)
C
      LOGICAL PACK2E
C
      DIMENSION WOCC3(NOC,NOC),PAI(NVIR,NOC),
     *          WOCC(NOC,NOC),WAI(NVIR,NOC),C(NBF,NBF)
      DIMENSION XLMN(NBF,NBF),PMN(NBF2),E(*),X(NBF,NBF),
     *          Y(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
C
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
C
C     ---- ADD FINAL TERMS TO THE ENERGY WEIGHTED DENSITY ----
C
      CALL VCLR(WOCC3,1,NOC*NOC)
C
C     ---- FORM FOCK-LIKE MATRIX ----
C
      CALL FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
C     ---- TRANSFORM THE FOCK LIKE MATRIX G(MU,NU) TO MO BASIS ----
C
      CALL MRARBR(Y,NBF,NBF,NBF,C,NBF,NOC,X,NBF)
      CALL DGEMM('T','N',NOC,NOC,NBF,ONE,C,NBF,X,NBF,
     *           ONE,WOCC3,NOC)
C
C     ---- COMPLETE THE ENERGY WEIGHTED DENSITY IN THE MO BASIS ----
C
      DO 50 J=1,NOC
         DO 60 I=1,NOC
            WOCC(I,J)=WOCC(I,J) - WOCC3(I,J)
   60    CONTINUE
   50 CONTINUE
C
C     ---- ADD ENERGY TERM TO WAI ----
C
      DO 125 I=1,NOC
         DO 120 J=1,NVIR
            WAI(J,I) = WAI(J,I) - E(I)*PAI(J,I)
  120    CONTINUE
  125 CONTINUE
C
      RETURN
      END
C*MODULE MP2GRD  *DECK AOCPHF
      SUBROUTINE AOCPHF(PMN,XX,Y,XLMN,BUF,IBUF,XLAI,PAI,PBJ,NOC,
     *                  NVIR,NBF,C,ERR,WRK,A,X,IPVT,B,IODIIS,
     *                  E,NROT,MAXIT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,NOTOPN
C
      DIMENSION PAI(NVIR,NOC), PBJ(NVIR,NOC), PMN(*)
      DIMENSION XX(NBF,NBF),   XLMN(*), XLAI(NVIR,NOC)
      DIMENSION Y(NBF,NBF),    BUF(NINTMX), IBUF(NINTMX)
      DIMENSION C(NBF,*), E(NBF), ERR(NVIR,NOC)
      DIMENSION WRK(NROT),A(MAXDII,MAXDII),X(MAXIT),IPVT(MAXIT),B(*)
      DIMENSION IODIIS(2*MAXDII)
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, CVGTOL=1.0D-10)
C
      ITDIIS = 1
      NOCP   = 0
      MAXIO  = 2*MAXDII
      MAXIT2 = (MAXIT*MAXIT+MAXIT)/2
      ERDIIS = ZERO
      NOTOPN = .TRUE.
C
C        SET UP AN INITIAL GUESS AT RESPONSE VECTORS
C
      CALL DCOPY(NROT,XLAI,1,PBJ,1)
      DO 60 NI=1,NOC
         DO 50 NA=1,NVIR
            PBJ(NA,NI)=PBJ(NA,NI)/(E(NA+NOC)-E(NI))
   50    CONTINUE
   60 CONTINUE
C
C  START DIIS ITERATIONS
C
      DELTA = ZERO
      ITER   = 0
  100 CONTINUE
      ITER=ITER+1
C
C     ----- FORM D(LAMDA,SIGMA) WHICH WE CALL PMN -----
C     ACCORDING TO EQUATION 17 OF THE FRISCH PAPER.
C
C        BEGIN TRANSFORMING -PBJ- TO AO BASIS AS -PMN-
C
      CALL VCLR(Y,1,NBF*NBF)
C
      CALL MRARTR(PBJ,NVIR,NVIR,NOC,C,NBF,NBF,XX,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,TWO,C(1,NOC+1),NBF,XX,NBF,
     *     ONE,Y,NBF)
C
      MUNU = 0
      DO MU=1,NBF
         DO NU=1,MU
            MUNU=MUNU+1
            PMN(MUNU) = Y(MU,NU) + Y(NU,MU)
         ENDDO
      ENDDO
C
C           NOW FORM FOCK-LIKE MATRIX
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
C     ---- TRANSFORM THE FOCK-LIKE MATRIX G(MU,NU) TO MO BASIS ----
C
      CALL VCLR(PAI,1,NROT)
C
      CALL MRARBR(Y,NBF,NBF,NBF,C,NBF,NOC,XX,NBF)
      CALL DGEMM('T','N',NVIR,NOC,NBF,ONE,C(1,NOC+1),NBF,XX,NBF,
     *           ONE,PAI,NVIR)
C
C     ---- CALCULATE PAI USING EQUATION 10 IN FRISCH'S PAPER ----
C
      DO 700 NI=1,NOC
         DO 710 NA=1,NVIR
            PAI(NA,NI)=(XLAI(NA,NI)-PAI(NA,NI))/(E(NA+NOC)-E(NI))
  710    CONTINUE
  700 CONTINUE
C
      IF(ITER.GE.2) THEN
C
C     --- DETERMINE ERROR VECTOR FOR DIIS ----
C
         DO 800 NI=1,NOC
            DO 810 NA=1,NVIR
               ERR(NA,NI) = PAI(NA,NI) - PBJ(NA,NI)
  810       CONTINUE
  800    CONTINUE
C
C     ---- CARRY OUT DIIS ----
C
         ETHSAV = ETHRSH
         ETHRSH = MAX(5.0D+00,ETHRSH)
         CALL CPDIIS(IW,ITDIIS,PAI,ERR,WRK,A,X,
     *               IPVT,B,IODIIS,NOCP,MAXIT,MAXIT2,
     *               MAXIO,ERDIIS,NOTOPN,NROT,NOC,NVIR,1)
         ETHRSH = ETHSAV
      END IF
C
C     ---- CHECK CONVERGENCE OF SOLUTION ----
C
      DELPRV = DELTA
      DELTA=ZERO
      DO 900 NI=1,NOC
         DO 910 NA=1,NVIR
            DELTA=DELTA + (PAI(NA,NI)-PBJ(NA,NI))**2
  910    CONTINUE
  900 CONTINUE
      DELTA=SQRT(DELTA/NROT)
      IF(NPRINT.NE.-5) WRITE(IW,9020) DELTA
C
C        NEXT IS CONVERGED EXIT
C
      IF(DELTA.LT.CVGTOL  .AND.  ITER.GE.2) THEN
         CALL RACLOS(IRAF,'DELETE')
         NOTOPN = .TRUE.
         WRITE(IW,9000) ITER
         IF(DBUG) THEN
            WRITE(IW,9030)
            CALL PRSQ(PAI,NOC,NVIR,NVIR)
         END IF
         RETURN
      END IF
C
      CALL DCOPY(NROT,PAI,1,PBJ,1)
C
C        CONTINUE WORKING ON CONVERGENCE
C
      IF(ITER.LT.MAXIT)  GO TO 100
C
C        THE CPHF HAS FAILED TO CONVERGE
C
      CALL RACLOS(IRAF,'DELETE')
      NOTOPN = .TRUE.
      WRITE(IW,9000) ITER
      WRITE(IW,9010) DELTA,DELPRV,CVGTOL
      CALL ABRT
      RETURN
C
 9000 FORMAT(1X,'.....',I4,' CPHF DIIS ITERATIONS CARRIED OUT .....'/)
 9010 FORMAT(1X,'AOCPHF: TOO MANY ITERATIONS IN THE CPHF'/
     *       1X,'   FINAL ITERATION Z-MATRIX ERROR=',1P,E13.6/
     *       1X,'PREVIOUS ITERATION Z-MATRIX ERROR=',1P,E13.6/
     *       1X,' CPHF CONVERGENCE CRITERION IS',1P,E13.6/
     *       1X,'PLEASE RECONSIDER:'/
     *       5X,'IF THE SCF HAS CONVERGED TO THE RIGHT ORBITAL',
     *          ' OCCUPANCY OR NOT'/
     *       5X,'IF THIS CHEMICAL SYSTEM IS A SINGLE',
     *          ' REFERENCE WAVEFUNCTION OR NOT.'//
     *       1X,'IF YOU BELIEVE A SINGLE REFERENCE WAVEFUNCTION IS',
     *          ' APPROPRIATE,'/
     *       1X,'AND IF AND ONLY IF THE CPHF ITERATIONS ABOVE SEEM TO',
     *          ' BE CONVERGING,'/
     *       1X,'INCREASE -MAXITC- (UNDOCUMENTED) IN $MP2 OR $CIS.')
 9020 FORMAT(1X,'MAXIMUM Z-VECTOR ERROR=',1P,E13.6)
 9030 FORMAT(/' PAI(VIRT,OCC)')
      END
C*MODULE SCFLIB  *DECK FLMAT
      SUBROUTINE FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E,SCHWRZ,DIRSCF,FDIFF,GOPARR,DSKWRK,MASWRK,TDSKWRK
C
      DIMENSION PMN(*), XLMN(*)
      DIMENSION Y(NBF,NBF), BUF(NINTMX), IBUF(NINTMX)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXAO=8192)
C
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (HALF=0.5D+00)
C
      DATA RHF/8HRHF     /
C
C           NOW FORM FOCK-LIKE MATRIX
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      NINT = 0
      NSCHWZ = 0
      L1 = NBF
      L2 = (NBF*NBF+NBF)/2
C
      IF(DIRSCF) THEN
         SCHWRZ=ISCHWZ.GT.0
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
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
         CALL VALFM(LOADFM)
         LGHOND = LOADFM + 1
         LXINTS = LGHOND + MAXG
         LDSH   = LXINTS + NSH2
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
         IF(SCHWRZ) THEN
            DUMMY = 0.0D+00
            CALL SHLDEN(RHF,PMN,DUMMY,DUMMY,XX(LDSH),IA,L1,L2,NSH2,1)
            CALL DAREAD(IDAF,IODA,XX(LXINTS),NSH2,54,0)
         END IF
         CALL VCLR(XLMN,1,L2)
         CALL TWOEI(RHF,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     *              INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *              DUMMY,DUMMY,IDUMM,1,XX(LXINTS),NSH2,XX(LGHOND),
     *              MAXG,XX(LDDIJ),IA,PMN,XLMN,DUMMY,DUMMY,XX(LDSH),
     *              DUMMY,DUMMY,1)
         CALL DSCAL(L2,HALF,XLMN,1)
         II=0
         DO 210 I=1,L1
            II = II+I
            XLMN(II) = XLMN(II) + XLMN(II)
  210    CONTINUE
C
         CALL RETFM(NEED)
      ELSE
         TDSKWRK = DSKWRK
         DSKWRK  = .TRUE.
         CALL SEQREW(IJK)
         CALL HSTAR(PMN,XLMN,BUF,IBUF,NINTMX,IA,NOPK,.FALSE.)
         DSKWRK = TDSKWRK
      END IF
C
C     --- SUM UP PARTS TO MAKE FOCK-LIKE MATRIX ---
C
      IF(GOPARR) THEN
         CALL DDI_GSUMI(2311,NINT,1)
         CALL DDI_GSUMI(2312,NSCHWZ,1)
         CALL DDI_GSUMF(2315,XLMN,L2)
      END IF
C
C     ----- SYMMETRIZE FOCK-LIKE MATRIX -----
C
      CALL SYMH(XLMN,Y,IA)
C     CALL PRTRI(XLMN,NBF)
C
      CALL VCLR(Y,1,NBF*NBF)
      MUNU = 0
      DO 20 MU=1,NBF
      DO 20 NU=1,MU
         MUNU=MUNU+1
         Y(MU,NU) = XLMN(MUNU)
   20    Y(NU,MU) = XLMN(MUNU)
C
      RETURN
      END
C        -         -         -         -         -         -         - C
C*MODULE MP2GRD  *DECK UMPGE2
      SUBROUTINE UMPGE2
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXRT=100, MXATM=2000, MXAO=8192)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,DIRSCF,FDIFF,DIRTRF,
     *        GOPARR,DSKWRK,MASWRK
C
      COMMON /DIRMEM/ IPQKL,IVEC,IDUM1,IDUM2,NUM,NOCC,NDIM,MINPQ,MAXPQ
      COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /ENRGMP/ EMP2,EMP3,EMP4
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NNUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWORD,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /MP2PRT/ SOME,OUT,DBUG
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
C     COSMO ADDITIONS FOR MP2 ANALYSIS
C
      LOGICAL ISEPS, USEPS
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /COSMO2/ QENUC,ELAST,EMP2COS,EMP2COS2,ETOTSMP,SAVESE,
     *                EMP2LAST,MP2TRIP,MP2ITER,MP2FACT
C
      DATA KRECA1,KRECA2,KRECA3,KRECA4,KRECA5,KRECA6
     *        /411,412,413,414,415,416/
      DATA KRECB1,KRECB2,KRECB3,KRECB4,KRECB5,KRECB6
     *        /421,422,423,424,425,426/
      DATA NFTIA,NFTIJ,NFT16/9,15,16/
      DATA NFTAB,NFTBA,NFTBB/41,42,43/
      DATA CHECK/8HCHECK   /, DBUGME,DEBUG/8HMP2GRAD ,8HDEBUG   /
C
C     ----- THIS ROUTINE CALCULATES THE ENERGY CONTRIBUTIONS TO THE
C     UMP2 GRADIENT.  AFTER THE ENERGY IS CALCULATED, SOME PART OF
C     DENSITY AND ENERGY WEIGHTED DENSITY ARE CALCULATED -----
C
C     ----- THIS CODE IS WRITTEN BY CHRISTINE AIKENS -----
C            IT IS MODELED AFTER THE MP2GE2 ROUTINE
C
      NOCA  = NOA - NACORE
      NOCB  = NOB - NBCORE
      NVIRA = NO  - NOA
      NVIRB = NO  - NOB
      NORBA = NO  - NACORE
      NORBB = NO  - NBCORE
      NORB  = NO
C
      NBF2  = (NBF*NBF+NBF)/2
      NBF3  = NBF*NBF
      NUM   = NBF
      NDIM = MAX(NBF2,NVIRA*NOA,NVIRB*NOB)
C
      SOME = NPRINT.NE.-5  .AND.  MASWRK
      OUT  = .FALSE.
      DBUG = MASWRK  .AND.  (EXETYP.EQ.DBUGME  .OR.  EXETYP.EQ.DEBUG)
      IF(MASWRK) WRITE(IW,9000)
      IF(SOME) WRITE(IW,9005) NACORE,NOCA,NBCORE,NOCB,NQMT,NBF
C
      CALL SEQOPN(NFTIA,'MOINTS','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFTAB,'AABB41','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFTBA,'BBAA42','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFTBB,'BBBB43','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFTIJ,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NFT16,'WORK16','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFTIA)
      CALL SEQREW(NFTAB)
      CALL SEQREW(NFTBA)
      CALL SEQREW(NFTBB)
C
C        ---- OPTION TO DEFEAT SYMMETRY USAGE ----
C
      IF(NOSYMT.EQ.1) CALL SYMOFF
C
      CALL GOTFM(NGOTMX)
      IF(NWORD.GT.0) NGOTMX=MIN(NWORD,NGOTMX)
C
C     LVECA:    NBF*NBF           - ALPHA ORBITAL COEFF (ACT + VIRT)
C     LVECB:    NBF*NBF           - BETA ORBITAL COEFF
C     IIRPA:    NBF               - ALPHA MOSIRP (IRREP LABELS FOR MO'S)
C     IIRPB:    NBF               - BETA MOSIRP
C     IWRK:     MAX(INITMX,NBF3)  - INTERMEDIATE DATA
C     IWRK2:    MAX(NINTMX,NBF3)  - INTERMEDIATE DATA
C     ITVEC:    NBF2              - OVERLAP INTEGRALS
C     IPOCC:    NBF3              - INTERMEDIATE COEFFS/Q MATRIX
C     ITMP:     NBF               - TEMPORARY DATA
C     LVEC2A:   NBF*NBF           - ALPHA ORBITAL COEFF
C     LVEC2B:   NBF*NBF           - BETA ORBITAL COEFF (CORE+ACT+VIRT)
C
      CALL VALFM(LOADFM)
      LVECA = 1     + LOADFM
      LVECB = LVECA + NBF3
      IIRPA = LVECB + NBF3
      IIRPB = IIRPA + NBF
      LAST  = IIRPB + NBF
      LVEC2A = LVECA
      LVEC2B = LVECB
      IF(NACORE.GT.0) THEN
        LVEC2A = LAST
        LAST  = LVEC2A + NBF3
      END IF
      IF(NBCORE.GT.0) THEN
        LVEC2B = LAST
        LAST = LVEC2B + NBF3
      END IF
      NEEDA = LAST  - LOADFM
      CALL GETFM(NEEDA)
C
C     ----- READ IN ORBITALS, AND ASSIGN SYMMETRY -----
C
      CALL VALFM(LOADFM)
      IWRK   = 1     + LOADFM
      IWRK2  = IWRK  + NBF
      ITVEC  = IWRK2 + NBF
      IPOCC  = ITVEC + NBF2
      ITMP   = IPOCC + NBF3
      LAST   = ITMP  + NBF
      NEEDB  = LAST  - LOADFM
      CALL GETFM(NEEDB)
      IF(EXETYP.EQ.CHECK) GO TO 90
C
      IF (NACORE.GT.0) THEN
C
C     ---- NEED ONLY ACTIVE ORBITAL COEFFICIENTS IN 1 INDEX ----
C     ---- TRANSFORMATION (INPENDING CALL TO MP2C1)         ----
C     ---- SHIFT OCCUPIED COEFFICIENTS TO LVECA AND LVECB   ----
C     ---- IF THERE ARE ALPHA AND BETA CORE ORBITALS        ----
C
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,15,0)
         DO 50 J=1,NORBA
            DO 60 I=1,NBF
               X(I+ NBF*(J-1)+LVECA-1)=X(I+NBF*(J-1+NACORE)+IPOCC-1)
  60        CONTINUE
  50     CONTINUE
C
C     ---- NEED ALL ORBITAL COEFFICIENTS TO ASSIGN SYMMETRY ----
C     ---- AND TRANSFORMATION OF REMAINING INDICES (UMPGC2) ----
C
         CALL DAREAD(IDAF,IODA,X(LVEC2A),NBF3,15,0)
         CALL DAREAD(IDAF,IODA,X(ITVEC),NBF2,12,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
C
         CALL TRFSYM(X(IWRK),X(IIRPA),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVEC2A),X(ITMP),IA,NQMT,NBF,NO,NBF)
C
      ELSE
C
         CALL DAREAD(IDAF,IODA,X(LVECA),NBF3,15,0)
         CALL DAREAD(IDAF,IODA,X(ITVEC),NBF2,12,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
C
         CALL TRFSYM(X(IWRK),X(IIRPA),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVECA),X(ITMP),IA,NQMT,NBF,NO,NBF)
      END IF
C
C     ---- DO THE SAME THING FOR BETA COEFFICIENTS ----
C
      IF (NBCORE.GT.0) THEN
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,19,0)
         DO 70 J=1,NORBB
            DO 80 I=1,NBF
               X(I+ NBF*(J-1)+LVECB-1)=X(I+NBF*(J-1+NBCORE)+IPOCC-1)
  80        CONTINUE
  70     CONTINUE
C
         CALL DAREAD(IDAF,IODA,X(LVEC2B),NBF3,19,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
C
         CALL TRFSYM(X(IWRK),X(IIRPB),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVEC2B),X(ITMP),IA,NQMT,NBF,NO,NBF)
C
      ELSE
C
         CALL DAREAD(IDAF,IODA,X(LVECB) ,NBF3,19,0)
         CALL DAREAD(IDAF,IODA,X(IPOCC),NBF3,45,0)
C
         CALL TRFSYM(X(IWRK),X(IIRPB),X(IWRK2),X(IPOCC),X(ITVEC),
     *               X(LVECB),X(ITMP),IA,NQMT,NBF,NO,NBF)
      END IF
C
   90 CONTINUE
      CALL RETFM(NEEDB)
C
C     -LVECA-   NBF*NBF            - ALPHA ORBITAL COEFF
C     -LVECB-   NBF*NBF            - BETA ORBITAL COEFF
C     -IIRPA-   NBF                - ALPHA MOSIRP
C     -IIRPB-   NBF                - BETA MOSIRP
C     THE ABOVE FOUR HAVE ALREADY BEEN ALLOCATED. NOW WE ALLOCATE FOR:
C     -IWRK -   MAX(NINTMX,NBF3,NBF*(NOCA+NOCB)) - INTERMEDIATE DATA
C     -IWRK2-   MAX(NINTMX,NBF3)                 - INTERMEDIATE DATA
C     -IWRK3-   MAX(NINTMX,NBF3,NBF*(NOCA+NOCB)) - INTERMEDIATE DATA
C     -IENGA-   NBF                - ALPHA ENERGIES (EA)
C     -IPOCCA-  NOA*NOA            - OCC-OCC BLOCK DENSITY (POCCA)
C     -IWOCCA-  NOA*NOA            - OCC-OCC E WEIGHTED DENS (WOCCA)
C     -IPVIRA-  NVIRA*NVIRA        - VIR-VIR DENSITY (PVIRA)
C     -IWVIRA-  NVIRA*NVIRA        - VIR-VIR E WEIGHTED DENS (WVIRA)
C     -ILAI1A-  NOA*NVIRA          - MP2 LAGRANGIAN LAI(I) FOR ALPHA
C     -ILAI2A-  NOA*NVIRA          - MP2 LAGRANGIAN LAI(II) FOR ALPHA
C     -ALSO HAVE ANALOGOUS BETA-
C     THE ABOVE ARE ALLOCATED IN NEEDC, AS IS REQUIRED DIRSCF MEMORY
C     THE FOLLOWING ARE ALLOCATED IN NEEDD:
C     -IPQKL -   NBF*NBF2*NOC      - (PQ|RJ) INTEGRALS
C     -IPQKL2-   NBF*NBF2*NOC      - (PQ|RJ) INTEGRALS
C     -IWRK4-   NOC*NVIR*NOC       - INTERMEDIATE DATA
C
      CALL VALFM(LOADFM)
      IWRK   = 1      + LOADFM
      IWRK2  = IWRK   + MAX(NINTMX,NBF3,NBF*(NOCA+NOCB))
      IWRK3  = IWRK2  + MAX(NINTMX,NBF3)
      IENGA  = IWRK3  + MAX(NINTMX,NBF3,NBF*(NOCA+NOCB))
      IENGB  = IENGA  + NBF
      IPOCCA = IENGB  + NBF
      IPOCCB = IPOCCA + NOA*NOA
      IWOCCA = IPOCCB + NOB*NOB
      IWOCCB = IWOCCA + NOA*NOA
      IPVIRA = IWOCCB + NOB*NOB
      IPVIRB = IPVIRA + NVIRA*NVIRA
      IWVIRA = IPVIRB + NVIRB*NVIRB
      IWVIRB = IWVIRA + NVIRA*NVIRA
      ILAI1A = IWVIRB + NVIRB*NVIRB
      ILAI1B = ILAI1A + NOA*NVIRA
      ILAI2A = ILAI1B + NOB*NVIRB
      ILAI2B = ILAI2A + NOA*NVIRA
      LAST   = ILAI2B + NOA*NVIRB
      LSTSAV = LAST
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
      NSH2 = (NSHELL*NSHELL+NSHELL)/2
      IF(DIRSCF) THEN
         IGHND  = LAST
         IXINTS = IGHND  + MAXG
         IDDIJ  = IXINTS + NSH2
         LAST   = IDDIJ  + 49*MXG2
      END IF
      NEEDC  = LAST - LOADFM
      CALL GETFM(NEEDC)
C
      IF(DIRSCF  .AND.  ISCHWZ.EQ.1)
     *     CALL DAREAD(IDAF,IODA,X(IXINTS),NSH2,54,0)
C
      NOCX  = MAX(NOCA,NOCB)
C
C     FIGURE OUT MEMORY REQUIREMENTS.  THE INCREMENT IS 2*(NBF2*NBF*2)
C     (PLUS THE FIRST TERM) TO ACCOUNT FOR THE TWO IPQKL MATRICES
C
      LEFT   = NGOTMX - NEEDA - NEEDC
      NMIN   = NEEDA + NEEDC + MAX(NVIRA*NOCA,NVIRB*NOCB) + NBF2*NBF*4
      NMAX   = NEEDA + NEEDC + MAX(NVIRA*NOCA*NOCA,NVIRB*NOCB*NOCB)
     *         +NBF2*NBF*NOCX*4
      NINCR  =                 MAX(NVIRA*NOCA,NVIRB*NOCB) + NBF2*NBF*4
      IF(SOME) WRITE(IW,9990) NMIN,NMAX,NINCR
C
      IF(NMIN.GT.NGOTMX) THEN
         IF(MASWRK) WRITE(IW,9997) NMIN-LEFT
         CALL ABRT
      END IF
C
      LPASS = MIN(NOCX,LEFT/NINCR)
      NPASS = (NOCX+LPASS-1)/LPASS
      LPASS = (NOCX+NPASS-1)/NPASS
C
      NO2AA  = 0
      NO2AB  = 0
      NO2BA  = 0
      NO2BB  = 0
C
      DO 102 I=IPOCCA,LSTSAV
  102 X(I)=ZERO
C
      CALL VALFM(LOADFM)
      IVEC   = 1      + LOADFM
      IPQKL  = IVEC   + NBF*LPASS*2
      IPQKL2 = IPQKL  + NBF*NBF2*LPASS*2
      IWRK4  = IPQKL2 + NBF*NBF2*LPASS*2
      LAST   = IWRK4  + LPASS*MAX(NVIRA*NOA,NVIRB*NOB)
      NEEDD  = LAST  - LOADFM
      CALL GETFM(NEEDD)
      NTOT  = NEEDA + NEEDC + NEEDD
      IF(SOME) WRITE(IW,9999) NGOTMX,NTOT,LPASS,NPASS,TOL
      IF(EXETYP.EQ.CHECK) GO TO 600
C
C        FOR THIS I2CASE=3 DIRECT TEI TRANSFORMATION, THE /DIRMEM/
C        COMMON MUST READ IN -IPQKL- AND -IVEC- (WHICH POINTS AT
C        THE FIRST MO OF EACH PARTICULAR PASS)
C
      IOCA=1
      IOCB=1
      IPASS=0
      CALL TSECND(TIM1)
C
C     SET UP A LOOP TO TRANSFORM ALL (PQ|RS) TO (AB|CJ) AND FORM
C     ONE ELECTRON DENSITY MATRICES
C
  100 CONTINUE
      IPASS=IPASS+1
      JOCA=MIN(NOCA,IOCA+LPASS-1)
      NOCCA=JOCA-IOCA+1
      JOCB=MIN(NOCB,IOCB+LPASS-1)
      NOCCB=JOCB-IOCB+1
      NOCC2=NOCCA+NOCCB
C
      CALL UMP2TR(X(IPQKL),X(IPQKL2),X(LVEC2A),X(LVEC2B),X(LVECA),
     *           X(LVECB),X(IVEC),X(IENGA),X(IENGB),X(IIRPA),X(IIRPB),
     *           X(IWRK4),
     *           X(IWRK),X(IWRK2),X(IWRK3),X(IWRK),X(IWRK3),X(IWRK2),
     *           NOCA,NVIRA,NOCB,NVIRB,NOCCA,NOCCB,NOCC2,NBF,NBF2,NDIM,
     *           NACORE,NBCORE,IOCA,IOCB,NORB,
     *           NOA,NOB,IJK,TOL,
     *           NO2AA,NFTIA,NO2AB,NFTAB,NO2BA,NFTBA,NO2BB,NFTBB,
     *           X(IGHND),X(IXINTS),X(IDDIJ),NINTMX,MAXG,NSH2,MXG2,
     *           X(IPOCCA),X(IWOCCA),X(IPVIRA),X(IWVIRA),X(ILAI1A),
     *           X(ILAI2A),X(IPOCCB),X(IWOCCB),X(IPVIRB),X(IWVIRB),
     *           X(ILAI1B),X(ILAI2B))
C
      IF(SOME) THEN
         CALL TSECND(TIM2)
         TIME = TIM2 - TIM1
         TIM1 = TIM2
         WRITE(IW,9980) IPASS,TIME
         CALL FLSHBF(IW)
      END IF
      IOCB = JOCB+1
      IOCA = JOCA+1
      IF(IPASS.LT.NPASS) GO TO 100
      IF(SOME) WRITE(IW,9985)
C
C     LOOP IS DONE.  WE NOW HAVE ALL (AB|CJ) AND ALPHA AND BETA MATRICES
C     GET UMP2 ENERGY AND ADD ENERGY CONTRIBUTIONS TO -WOCC- AND -WVIR-
C
      CALL UMPGEN(E2,X(IENGA),X(IENGB),X(IPOCCA),X(IPOCCB),X(IWOCCA),
     *            X(IWOCCB),X(IPVIRA),X(IPVIRB),X(IWVIRA),X(IWVIRB),
     *            NOA,NOB,NVIRA,NVIRB)
C
      EHF  = ESCF
      EMP2 = EHF+E2
C
C     COSMO CHANGES
C
      IF(ISEPS) THEN
         WRITE(IW,*)"INSIDE UMPGE2:"
         WRITE(IW,*)"VALUE OF ISEPS IN MP2:",ISEPS
         IF(MP2TRIP.EQ.0) EMP2COS=EMP2
         MP2TRIP=1
         WRITE(IW,7802) EHF,E2,EMP2
      ELSE
         WRITE(IW,7800) EHF,E2,EMP2
      ENDIF
C
C
      IF(DBUG) THEN
         WRITE(IW,8999)
         CALL PRSQ(X(IPOCCA),NOA,NOA,NOA)
         WRITE(IW,8998)
         CALL PRSQ(X(IPVIRA),NVIRA,NVIRA,NVIRA)
         WRITE(IW,8997)
         CALL PRSQ(X(IWOCCA),NOA,NOA,NOA)
         WRITE(IW,8996)
         CALL PRSQ(X(IWVIRA),NVIRA,NVIRA,NVIRA)
         WRITE(IW,8995)
         CALL PRSQ(X(ILAI1A),NOA,NVIRA,NVIRA)
         WRITE(IW,8994)
         CALL PRSQ(X(ILAI2A),NOA,NVIRA,NVIRA)
         WRITE(IW,8989)
         CALL PRSQ(X(IPOCCB),NOB,NOB,NOB)
         WRITE(IW,8988)
         CALL PRSQ(X(IPVIRB),NVIRB,NVIRB,NVIRB)
         WRITE(IW,8987)
         CALL PRSQ(X(IWOCCB),NOB,NOB,NOB)
         WRITE(IW,8986)
         CALL PRSQ(X(IWVIRB),NVIRB,NVIRB,NVIRB)
         WRITE(IW,8985)
         CALL PRSQ(X(ILAI1B),NOB,NVIRB,NVIRB)
         WRITE(IW,8984)
         CALL PRSQ(X(ILAI2B),NOB,NVIRB,NVIRB)
      END IF
C
      IF(SOME) WRITE(IW,9996) NO2AA,NFTIA,NO2AB,NFTAB,
     *                        NO2BA,NFTBA,NO2BB,NFTBB
C
C     ----- SAVE ON FILE -DA10- -----
C
C     -POCCA- ON RECORD -KRECA1-
C     -PVIRA- ON RECORD -KRECA2-
C     -WAIA - ON RECORD -KRECA3-
C     -LAIA - ON RECORD -KRECA4-
C     -WOCCA- ON RECORD -KRECA5-
C     -WVIRA- ON RECORD -KRECA6-
C     -POCCB- ON RECORD -KRECB1-
C     -PVIRB- ON RECORD -KRECB2-
C     -WAIB - ON RECORD -KRECB3-
C     -LAIB - ON RECORD -KRECB4-
C     -WOCCB- ON RECORD -KRECB5-
C     -WVIRB- ON RECORD -KRECB6-
C         IN OTHER PLACES, WE WILL SAVE OTHER RECORDS:
C     -P(2)A- ON RECORD -KRECA7-
C     -PHFA-  ON RECORD -KRECA8-
C     -W(2)A- ON RECORD -KRECA9-
C     -P(2)B- ON RECORD -KRECB7-
C     -PHFB-  ON RECORD -KRECB8-
C     -W(2)B- ON RECORD -KRECB9-
C
  600 CONTINUE
C
      CALL DAWRIT(IDAF,IODA,X(IPOCCA),NOA  *NOA  ,KRECA1,0)
      CALL DAWRIT(IDAF,IODA,X(IPVIRA),NVIRA*NVIRA,KRECA2,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI2A),NVIRA*NOA  ,KRECA3,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI1A),NVIRA*NOA  ,KRECA4,0)
      CALL DAWRIT(IDAF,IODA,X(IWOCCA),NOA  *NOA  ,KRECA5,0)
      CALL DAWRIT(IDAF,IODA,X(IWVIRA),NVIRA*NVIRA,KRECA6,0)
      CALL DAWRIT(IDAF,IODA,X(IPOCCB),NOB  *NOB  ,KRECB1,0)
      CALL DAWRIT(IDAF,IODA,X(IPVIRB),NVIRB*NVIRB,KRECB2,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI2B),NVIRB*NOB  ,KRECB3,0)
      CALL DAWRIT(IDAF,IODA,X(ILAI1B),NVIRB*NOB  ,KRECB4,0)
      CALL DAWRIT(IDAF,IODA,X(IWOCCB),NOB  *NOB  ,KRECB5,0)
      CALL DAWRIT(IDAF,IODA,X(IWVIRB),NVIRB*NVIRB,KRECB6,0)
C
      IF(NOSYMT.EQ.1) CALL SYMON
C
C     ----- ENERGY IS CALCULATED, RETURN TO THE CALLING PROGRAM -----
C
      CALL RETFM(NEEDD)
      CALL RETFM(NEEDC)
      CALL RETFM(NEEDA)
      RETURN
C
 8999 FORMAT(/' ----- P2(OCC,OCC) (ALPHA,ALPHA)    -----'/)
 8998 FORMAT(/' ----- P2(VIR,VIR) (ALPHA,ALPHA)    -----'/)
 8997 FORMAT(/' ----- W2(OCC,OCC) (ALPHA,ALPHA)    -----'/)
 8996 FORMAT(/' ----- W2(VIR,VIR) (ALPHA,ALPHA)    -----'/)
 8995 FORMAT(/' ----- LAGRANGIAN (I)  (ALPHA,ALPHA)-----'/)
 8994 FORMAT(/' ----- W2(VIR,OCC) (I) (ALPHA,ALPHA)-----'/)
 8989 FORMAT(/' ----- P2(OCC,OCC)   (BETA,BETA)    -----'/)
 8988 FORMAT(/' ----- P2(VIR,VIR)   (BETA,BETA)    -----'/)
 8987 FORMAT(/' ----- W2(OCC,OCC)   (BETA,BETA)    -----'/)
 8986 FORMAT(/' ----- W2(VIR,VIR)   (BETA,BETA)    -----'/)
 8985 FORMAT(/' ----- LAGRANGIAN (I)   (BETA,BETA) -----'/)
 8984 FORMAT(/' ----- W2(VIR,OCC) (I)  (BETA,BETA) -----'/)
 9000 FORMAT(/1X,16('-'),3X,29(1H-)/
     *        1X,'UHF-MP2 GRADIENT',3X,
     *           'PROGRAM WRITTEN BY C.M.AIKENS'/
     *        1X,16('-'),3X,29(1H-))
 9005 FORMAT(/1X,'# ALPHA CORE ORBITALS        =',I5/
     *        1X,'# ALPHA OCCUPIED ORBITALS    =',I5/
     *        1X,'# BETA CORE ORBITALS         =',I5/
     *        1X,'# BETA OCCUPIED ORBITALS     =',I5/
     *        1X,'# A OR B MOLECULAR ORBITALS  =',I5/
     *        1X,'# BASIS FUNCTIONS            =',I5)
 7800 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/)
 7802 FORMAT(/1X,'RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE'/
     *       13X,'  ECOSMO(0)               =',1X,F20.10/
     *       13X,'  E(1)                    =',9X,'0.0'/
     *       13X,'  E(2)                    =',1X,F20.10/
     *       13X,'E(COSMO-MP2)              =',1X,F20.10/)
C
 9999 FORMAT(/' UMP2 ENERGY/GRADIENT TEI TRANSFORMATION'/
     *        ' # OF WORDS AVAILABLE                = ',I10/
     *        ' # OF WORDS USED                     = ',I10/
     *        ' # OF ALPHA AND BETA ORBITALS/PASS   = ',I10/
     *        ' # OF INTEGRAL PASSES                = ',I10/
     *        ' THRESHOLD FOR KEEPING TRANSFORMED 2E- INTEGRALS = ',
     *        1P,E10.3/)
 9997 FORMAT(/' NOT ENOUGH MEMORY IN -UMP2- TRANSFORMATION'/
     *        ' INCREASE MEMORY SIZE BY AT LEAST ',I10,' WORDS')
 9996 FORMAT(1X,'DONE WITH -UMP2- GRADIENT TEI TRANFORMATION'/
     *     1X,I10,' NON-ZERO AAAA [AI|BJ] INTEGRALS WRITTEN TO FILE',I3/
     *     1X,I10,' NON-ZERO AABB [AI|BJ] INTEGRALS WRITTEN TO FILE',I3/
     *     1X,I10,' NON-ZERO BBAA [AI|BJ] INTEGRALS WRITTEN TO FILE',I3/
     *     1X,I10,' NON-ZERO BBBB [AI|BJ] INTEGRALS WRITTEN TO FILE',I3)
 9980 FORMAT(1X,'PASS #',I4,' TOOK',F10.2,' SECONDS.')
 9985 FORMAT(1X,'DONE WITH PARTIAL TRANSFORMATION (PQ|RS) TO (IA|JB)')
 9990 FORMAT(/1X,'MEMORY REQUIREMENTS FOR UMP2 TRANSFORMATION'/
     *       1X,' MINIMUM=',I10,' WORDS, USING 1 ALPHA AND BETA',
     *          ' ORBITAL PER PASS'/
     *       1X,' MAXIMUM=',I10,' WORDS, USING ONLY 1 AO INTEGRAL PASS'/
     *       1X,'EACH ADDITIONAL',I10,' WORDS REDUCES PASS COUNT BY 1.')
      END
C
C*MODULE MP2GRD  *DECK UMP2TR
      SUBROUTINE UMP2TR(XPQRJ1,XPQRJ2,VECA,VECB,CA,CB,C,EA,EB,MOSIRA,
     *                  MOSIRB,T2,X,Y,Z,XX,ZZ,IX,NOCA,NVIRA,NOCB,NVIRB,
     *                  NOCCA,NOCCB,NOCC2,NBF,NBF2,NDIM,
     *                  NACORE,NBCORE,IOCA,IOCB,NORB,NOA,
     *                  NOB,NFT,TOL,
     *                  NO2AA,NFTIA,NO2AB,NFTAB,NO2BA,NFTBA,NO2BB,NFTBB,
     *                  GHONDO,XINTS,DDIJ,NINTMX,MAXG,NSH2,MXG2,
     *                  POCCA,WOCCA,PVIRA,WVIRA,XLAI1A,XLAI2A,
     *                  POCCB,WOCCB,PVIRB,WVIRB,XLAI1B,XLAI2B)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION XPQRJ1(*),XPQRJ2(*),VECA(NBF,*),VECB(NBF,*),
     *          CA(NBF,*),CB(NBF,*),EA(*),EB(*),MOSIRA(*),MOSIRB(*),
     *          T2(*),C(NBF,*),X(NBF,*),Y(NBF,*),Z(NBF,*),XX(*),ZZ(*),
     *          IX(NINTMX),GHONDO(MAXG),XINTS(NSH2),DDIJ(49,MXG2)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      LOGICAL DIRSCF,ABEL,ABELPT,FDIFF
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
C     ----- THIS ROUTINE TRANSFORMS INTEGRALS, THEN FILLS -----
C     ----- PART OF DENSITY MATRICES                         -----
C
      ABEL=ABELPT()
      NTSAVE = NT
      IF(DIRSCF  .AND.  .NOT.ABEL) NT=1
C
C     ------ COPY -ALPHA- AND -BETA- ORBITALS INTO X MATRIX -----
C
      DO 10 I =1,NOCCA
      DO 10 MP=1,NBF
   10    C(MP,I      ) = CA(MP,IOCA+I-1)
      DO 20 I =1,NOCCB
      DO 20 MP=1,NBF
   20    C(MP,I+NOCCA) = CB(MP,IOCB+I-1)
C
C     ----- (PQ|RS) -> (PQ|RJ) -ALPHA- AND -BETA- -----
C
      CALL MP2C1(NFT,XPQRJ1,XX,IX,C,NBF,NOCC2,NDIM,
     *           GHONDO,XINTS,DDIJ,MAXG,NSH2)
C
      CALL UMPGC2(XPQRJ1,XPQRJ2,VECA,VECB,X,Y,Z,XX,ZZ,MOSIRA,MOSIRB,
     *            NBF,NBF2,NACORE,NBCORE,NOCCA,NOCCB,IOCA,IOCB,
     *            NOCC2,NOCA,NOCB,NORB,TOL,
     *            NO2AA,NFTIA,NO2AB,NFTAB,NO2BA,NFTBA,NO2BB,NFTBB)
C
      NT=NTSAVE
C
C    ----- READ IN ALPHA AND BETA ENERGIES -----
C
      CALL DAREAD(IDAF,IODA,EA,NBF,17,0)
      CALL DAREAD(IDAF,IODA,EB,NBF,21,0)
C
C    ----- CREATE P(AA),W(AA),L(AA) MATRICES -----
C
      CALL UP2W2N(XPQRJ1,T2,TOL,EA,EB,POCCA,WOCCA,PVIRA,WVIRA,
     *            XLAI1A,XLAI2A,NOA,NOB,NVIRA,NVIRB,NOCCA,NOCCB,IOCA,
     *            IOCB,NBF2,NBF,NACORE,NBCORE,NOCA,NOCC2,1)
C
C    ----- CREATE P(BB),W(BB),L(BB) MATRICES -----
C
      CALL UP2W2N(XPQRJ2,T2,TOL,EB,EA,POCCB,WOCCB,PVIRB,WVIRB,
     *            XLAI1B,XLAI2B,NOB,NOA,NVIRB,NVIRA,NOCCB,NOCCA,IOCB,
     *            IOCA,NBF2,NBF,NBCORE,NACORE,NOCB,NOCC2,2)
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UMPGC2
      SUBROUTINE UMPGC2(XPQRJ1,XPQRJ2,CA,CB,X,Y,Z,XX,ZZ,MOSIRA,MOSIRB,
     *                  NBF,NBF2,NACORE,NBCORE,NOCCA,NOCCB,IOCA,IOCB,
     *                  NOCC,NOCA,NOCB,NORB,TOL,
     *                  NO2AA,NFTIA,NO2AB,NFTAB,NO2BA,NFTBA,NO2BB,NFTBB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      LOGICAL SOME,OUT,DBUG,SYM2EI
C
      COMMON /MP2PRT/ SOME,OUT,DBUG
C
      DIMENSION XPQRJ1(NOCC,NBF2,NBF),XPQRJ2(NOCC,NBF2,NBF),
     *          CA(NBF,NORB),CB(NBF,NORB),X(NBF,*),Y(NBF,*),Z(NBF,*),
     *          XX(*),ZZ(*),MOSIRA(NORB),MOSIRB(NORB)
C
C     IN THIS ROUTINE, NORB IS THE TOTAL # OF ORBITALS
C     (CORE + ACTIVE + VIRTUAL)
C
      NORB2 = (NORB*NORB+NORB)/2
C
C     SECOND TRANSFORMATION GOES ONLY OVER VIRTUAL ORBS
C     (SET UP VARIABLES FOR THIS- NTRANA IS JUST NVIRA)
C
      NPA1   = NACORE+NOCA+1
      NTRANA = NORB-NACORE-NOCA
      NPB1   = NBCORE+NOCB+1
      NTRANB = NORB-NBCORE-NOCB
C
C     FINAL THREE INDICES OF THE INTEGRAL TRANSFORMATION
C     ROUTINE MODELED AFTER MP2GC2
C
C     ----- (PQ|RJ) TO (PQ|AJ) -----
C
C     FIRST DO ALPHA ALPHA COMBINATION
C
      DO 20 MPQ=1,NBF2
         DO 15 MR=1,NBF
            DO 10 MJ=1,NOCCA
               X(MR,MJ) = XPQRJ1(MJ,MPQ,MR)
   10       CONTINUE
   15    CONTINUE
C
         CALL MRTRBR(CA(1,NPA1),NBF,NBF,NTRANA,X,NBF,NOCCA,Y(NPA1,1),
     *               NBF)
C
         DO 30 MA=NPA1,NORB
            DO 25 MJ=1,NOCCA
               XPQRJ2(MJ,MPQ,MA) = Y(MA,MJ)
   25       CONTINUE
   30    CONTINUE
C
C     NEXT DO BETA BETA COMBINATION
C
         DO 45 MR=1,NBF
            DO 40 MJ=1,NOCCB
               X(MR,MJ) = XPQRJ1(MJ+NOCCA,MPQ,MR)
   40       CONTINUE
   45    CONTINUE
C
         CALL MRTRBR(CB(1,NPB1),NBF,NBF,NTRANB,X,NBF,NOCCB,Y(NPB1,1),
     *               NBF)
C
         DO 55 MA=NPB1,NORB
            DO 50 MJ=1,NOCCB
               XPQRJ2(MJ+NOCCA,MPQ,MA) = Y(MA,MJ)
   50       CONTINUE
   55    CONTINUE
   20 CONTINUE
C
C     ----- (PQ|AJ) TO (BC|AJ) FOR A J = ALPHA ALPHA-----
C
      DO 80 MJ=1,NOCCA
      DO 75 MA=NPA1,NORB
         JKT1 = NFTIA
         JKT2 = NFTAB
         MPQ = 0
         DO 90 MP=1,NBF
            DO 85 MQ=1,MP
               MPQ = MPQ+1
               X(MP,MQ) = XPQRJ2(MJ,MPQ,MA)
               X(MQ,MP) = X(MP,MQ)
   85       CONTINUE
   90    CONTINUE
         CALL MRARBR(X,NBF,NBF,NBF,CA,NBF,NORB,Y,NBF)
         CALL MRARBR(X,NBF,NBF,NBF,CB,NBF,NORB,Z,NBF)
         MPQ=0
         DO 120 MP=1,NORB
            DO 115 MQ=1,MP
               MPQ = MPQ+1
               IF(SYM2EI(MOSIRA(MP),MOSIRA(MQ),
     *                   MOSIRA(MA),MOSIRA(MJ+NACORE+IOCA-1))) THEN
                  XX(MPQ) = DDOT(NBF,Y(1,MP),1,CA(1,MQ),1)
                  IF(ABS(XX(MPQ)).LT.TOL) XX(MPQ) = ZERO
               ELSE
                  XX(MPQ) = ZERO
               END IF
               IF(SYM2EI(MOSIRB(MP),MOSIRB(MQ),
     *                   MOSIRA(MA),MOSIRA(MJ+NACORE+IOCA-1))) THEN
                  ZZ(MPQ) = DDOT(NBF,Z(1,MP),1,CB(1,MQ),1)
                  IF(ABS(ZZ(MPQ)).LT.TOL) ZZ(MPQ) = ZERO
               ELSE
                  ZZ(MPQ) = ZERO
               END IF
  115       CONTINUE
  120    CONTINUE
C
         DO 100 MPQ=1,NORB2
            XPQRJ1(MJ,MPQ,MA) = XX(MPQ)
            XPQRJ2(MJ,MPQ,MA) = ZZ(MPQ)
  100    CONTINUE
C
C        ----- STORE INTEGRALS ON DISK THAT WE WILL NEED AGAIN -----
C        ----- (WON'T NEED (VV|VO) AGAIN)                      -----
C
         MPQ=0
         DO 110 MP=1,NORB
            DO 105 MQ=1,MP
               MPQ=MPQ+1
               IF((MP.GT.NOCA+NACORE).AND.(MQ.GT.NOCA+NACORE))
     *            XX(MPQ) = ZERO
               IF(MP.GT.(NOCB+NBCORE).AND.MQ.GT.(NOCB+NBCORE))
     *            ZZ(MPQ) = ZERO
  105       CONTINUE
  110    CONTINUE
C
         CALL TRFWTM(JKT1,XX,NX,Y,NIX,NORB2,1,TOL)
         IF(DBUG.AND.NX.GT.0) THEN
            CALL TRFOUT(XX,NX,Y,NIX,NORB,MA,MJ+IOCA-1+NACORE)
         END IF
         CALL TRFWTM(JKT2,ZZ,NZ,Y,NIX,NORB2,1,TOL)
         IF(DBUG.AND.NZ.GT.0) THEN
            CALL TRFOUT(ZZ,NZ,Y,NIX,NORB,MA,MJ+IOCA-1+NACORE)
         END IF
         NO2AA = NO2AA+NX
         NO2AB = NO2AB+NZ
C
   75 CONTINUE
   80 CONTINUE
C
C     ----- (PQ|AJ) TO (BC|AJ) FOR A J = BETA BETA-----
C
      DO 180 MJ=1,NOCCB
      DO 175 MA=NPB1,NORB
         JKT1 = NFTBA
         JKT2 = NFTBB
         MPQ = 0
         DO 140 MP=1,NBF
            DO 135 MQ=1,MP
               MPQ = MPQ+1
               X(MP,MQ) = XPQRJ2(MJ+NOCCA,MPQ,MA)
               X(MQ,MP) = X(MP,MQ)
  135        CONTINUE
  140     CONTINUE
         CALL MRARBR(X,NBF,NBF,NBF,CA,NBF,NORB,Z,NBF)
         CALL MRARBR(X,NBF,NBF,NBF,CB,NBF,NORB,Y,NBF)
         MPQ=0
         DO 150 MP=1,NORB
            DO 145 MQ=1,MP
               MPQ = MPQ+1
               IF(SYM2EI(MOSIRB(MP),MOSIRB(MQ),
     *                   MOSIRB(MA),MOSIRB(MJ+NBCORE+IOCB-1))) THEN
                  XX(MPQ) = DDOT(NBF,Y(1,MP),1,CB(1,MQ),1)
                  IF(ABS(XX(MPQ)).LT.TOL) XX(MPQ) = ZERO
               ELSE
                  XX(MPQ) = ZERO
               END IF
               IF(SYM2EI(MOSIRA(MP),MOSIRA(MQ),
     *                   MOSIRB(MA),MOSIRB(MJ+NBCORE+IOCB-1))) THEN
                  ZZ(MPQ) = DDOT(NBF,Z(1,MP),1,CA(1,MQ),1)
                  IF(ABS(ZZ(MPQ)).LT.TOL) ZZ(MPQ) = ZERO
               ELSE
                  ZZ(MPQ) = ZERO
               END IF
  145        CONTINUE
  150     CONTINUE
C
         DO 160 MPQ=1,NORB2
            XPQRJ2(MJ+NOCCA,MPQ,MA) = XX(MPQ)
            XPQRJ1(MJ+NOCCA,MPQ,MA) = ZZ(MPQ)
  160     CONTINUE
C
C        ----- STORE INTEGRALS ON DISK THAT WE WILL NEED AGAIN -----
C        ----- (WON'T NEED (VV|VO) AGAIN)                      -----
C
         MPQ=0
         DO 170 MP=1,NORB
            DO 165 MQ=1,MP
               MPQ=MPQ+1
               IF(MP.GT.NOCB+NBCORE.AND.MQ.GT.NOCB+NBCORE)
     *            XX(MPQ) = ZERO
               IF(MP.GT.NOCA+NACORE.AND.MQ.GT.NOCA+NACORE)
     *            ZZ(MPQ) = ZERO
  165        CONTINUE
  170     CONTINUE
C
         CALL TRFWTM(JKT2,XX,NX,Y,NIX,NORB2,1,TOL)
         IF(DBUG.AND.NX.GT.0) THEN
            CALL TRFOUT(XX,NX,Y,NIX,NORB,MA,MJ+IOCB-1+NBCORE)
         END IF
         CALL TRFWTM(JKT1,ZZ,NZ,Y,NIX,NORB2,1,TOL)
         IF(DBUG.AND.NZ.GT.0) THEN
            CALL TRFOUT(ZZ,NZ,Y,NIX,NORB,MA,MJ+IOCB-1+NBCORE)
         END IF
         NO2BB = NO2BB+NX
         NO2BA = NO2BA+NZ
C
  175  CONTINUE
  180  CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UP2W2N
      SUBROUTINE UP2W2N(XPQRJ,T2,TOL,EA,EB,POCC,WOCC,PVIR,WVIR,
     *                 XLAI1,XLAI2,NOA,NOB,NVIRA,NVIRB,NOCCA,NOCCB,IOCA,
     *                 IOCB,NBF2,NBF,NACORE,NBCORE,NOCA,NOCC,IPASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
      PARAMETER (ZERO=0.0D+00, PT1=0.1D+00)
C
      COMMON /IJPAIR/ IA(MXAO)
C
      DIMENSION XPQRJ(NOCC,NBF2,NBF),T2(NVIRA,NOA,*)
      DIMENSION EA(*),EB(*),POCC(NOA,NOA),WOCC(NOA,NOA),
     *          PVIR(NVIRA,NVIRA),WVIR(NVIRA,NVIRA),
     *          XLAI1(NVIRA,NOA),XLAI2(NVIRA,NOA)
C
      TOL1 = TOL*PT1
C
C     THIS NEXT BIT ADJUSTS MEMORY STORAGE FOR XPQRJ
C
      IF (IPASS.EQ.1) THEN
         ISHIFT=0
         JSHIFT=NOCCA
      ELSE
         ISHIFT=NOCCB
         JSHIFT=0
      END IF
C
C    THIS ROUTINE IS BASED ON P2W2N FOR RHF-MP2 GRADIENTS
C
      DO 150 MB=1,NVIRA
C
C     ----- FIRST THE AMPLITUDES FROM A A A A SPIN CASE ------
C
      DO 104 MA=1,NVIRA
         MAJ = IA(MA+NOA)
         MBJ = IA(MB+NOA)
         DO 102 MJ=1,NOA
            MAJ = MAJ+1
            MBJ = MBJ+1
            DO 100 MK=1,NOCCA
               DJKAB = EA(MJ) + EA(MK+NACORE+IOCA-1) -
     *                 EA(MA+NOA) - EA(MB+NOA)
               TERM = XPQRJ(MK+ISHIFT,MAJ,MB+NOA) -
     *                XPQRJ(MK+ISHIFT,MBJ,MA+NOA)
               TERM = TERM/DJKAB
               IF(ABS(TERM).LT.TOL1) TERM = ZERO
               T2(MA,MJ,MK) = TERM
  100       CONTINUE
  102    CONTINUE
  104 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO POCC, WOCC, AND LAI(II) (FROM A A A A)
C
      DO 120 MJ=NACORE+1,NOA
         DO 118 MK=1,NOCCA
            DO 114 MA=1,NVIRA
C
               TERM = T2(MA,MJ,MK)
               IF(TERM.EQ.ZERO) GO TO 114
C
               MAI = IA(MA+NOA)
               DO 110 MI=1,NOA
                  MAI = MAI+1
                  MIJ = IA(MAX(MI,MJ)) + MIN(MI,MJ)
                  DIKAB = EA(MI) + EA(MK+NACORE+IOCA-1) - EA(MA+NOA) -
     *                    EA(MB+NOA)
                  XLAI2(MA,MI) = XLAI2(MA,MI) -
     *                           TERM*XPQRJ(MK+ISHIFT,MIJ,MB+NOA)
                  WOCC(MI,MJ) = WOCC(MI,MJ) -
     *                          TERM*XPQRJ(MK+ISHIFT,MAI,MB+NOA)
                  IF(MI.LE.NACORE) THEN
                     DIJ = EA(MJ) - EA(MI)
                     POCC(MI,MJ) = POCC(MI,MJ) +
     *                           TERM*XPQRJ(MK+ISHIFT,MAI,MB+NOA)/DIJ
                  ELSE
                     POCC(MI,MJ) = POCC(MI,MJ) -
     *                          TERM*XPQRJ(MK+ISHIFT,MAI,MB+NOA)/DIKAB
                  END IF
  110          CONTINUE
  114       CONTINUE
  118    CONTINUE
  120 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO PVIR, WVIR, AND LAI(I) (FROM A A A A)
C
      DO 140 MI=1,NOCA
         DO 138 MK=1,NOCCA
            DO 134 MC=1,NVIRA
C
               TERM = T2(MC,MI+NACORE,MK)
               IF(TERM.EQ.ZERO) GO TO 134
C
               DO 130 MA=1,NVIRA
                  MAI = IA(MA+NOA) + MI+NACORE
                  MAC = IA(MAX(MA,MC)+NOA) + MIN(MA,MC) + NOA
                  DIKAB = EA(MI+NACORE) + EA(MK+NACORE+IOCA-1) -
     *                    EA(MA+NOA) - EA(MB+NOA)
                  PVIR(MA,MC) = PVIR(MA,MC) +
     *                      TERM*XPQRJ(MK+ISHIFT,MAI,MB+NOA)/DIKAB
                  WVIR(MA,MC) = WVIR(MA,MC) -
     *                      TERM*XPQRJ(MK+ISHIFT,MAI,MB+NOA)
                  XLAI1(MA,MI+NACORE)= XLAI1(MA,MI+NACORE) +
     *                      TERM*XPQRJ(MK+ISHIFT,MAC,MB+NOA)
  130          CONTINUE
  134       CONTINUE
  138    CONTINUE
  140 CONTINUE
C
  150 CONTINUE
C
C     ----- NOW THE P AND W MATRICES HAVE ALL CONTRIBUTIONS  -----
C     ----- FROM AAAA.  NEXT WE DO AABB.                     -----
C
      DO 350 MB=1,NVIRB
C
C     ----- FIRST THE "AMPLITUDES" FROM A A B B SPIN CASE ------
C     THESE HAVE FORM (AJ|BK)/(EJ+EK-EA-EB), WHERE K AND B ARE BETA ORBS
C
      DO 304 MA=1,NVIRA
         MAJ = IA(MA+NOA)
         DO 302 MJ=1,NOA
            MAJ = MAJ+1
            DO 300 MK=1,NOCCB
               DJKAB = EA(MJ) + EB(MK+NBCORE+IOCB-1) -
     *                 EA(MA+NOA) - EB(MB+NOB)
               TERM = XPQRJ(MK+JSHIFT,MAJ,MB+NOB)
               TERM = TERM/DJKAB
               IF(ABS(TERM).LT.TOL1) TERM = ZERO
               T2(MA,MJ,MK) = TERM
  300       CONTINUE
  302    CONTINUE
  304 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO POCC, WOCC, AND LAI(II) (FROM A A B B)
C
      DO 320 MJ=NACORE+1,NOA
         DO 318 MK=1,NOCCB
            DO 314 MA=1,NVIRA
C
               TERM = T2(MA,MJ,MK)
               IF(TERM.EQ.ZERO) GO TO 314
C
               MAI = IA(MA+NOA)
               DO 310 MI=1,NOA
                  MAI = MAI+1
                  MIJ = IA(MAX(MI,MJ)) + MIN(MI,MJ)
                  DIKAB = EA(MI) + EB(MK+NBCORE+IOCB-1) - EA(MA+NOA) -
     *                    EB(MB+NOB)
                  XLAI2(MA,MI) = XLAI2(MA,MI) -
     *                           TERM*XPQRJ(MK+JSHIFT,MIJ,MB+NOB)
                  WOCC(MI,MJ) = WOCC(MI,MJ) -
     *                          TERM*XPQRJ(MK+JSHIFT,MAI,MB+NOB)
                  IF(MI.LE.NACORE) THEN
                     DIJ = EA(MJ) - EA(MI)
                     POCC(MI,MJ) = POCC(MI,MJ) +
     *                             TERM*XPQRJ(MK+JSHIFT,MAI,MB+NOB)/DIJ
                  ELSE
                     POCC(MI,MJ) = POCC(MI,MJ) -
     *                           TERM*XPQRJ(MK+JSHIFT,MAI,MB+NOB)/DIKAB
                  END IF
  310          CONTINUE
  314       CONTINUE
  318    CONTINUE
  320 CONTINUE
C
C     ----- NOW CONTRIBUTIONS TO PVIR, WVIR, AND LAI(I) (FROM A A B B)
C
      DO 340 MI=1,NOCA
         DO 338 MK=1,NOCCB
            DO 334 MC=1,NVIRA
C
               TERM = T2(MC,MI+NACORE,MK)
               IF(TERM.EQ.ZERO) GO TO 334
C
               DO 330 MA=1,NVIRA
                  MAI = IA(MA+NOA) + MI+NACORE
                  MAC = IA(MAX(MA,MC)+NOA) + MIN(MA,MC) + NOA
                  DIKAB = EA(MI+NACORE) + EB(MK+NBCORE+IOCB-1) -
     *                    EA(MA+NOA) - EB(MB+NOB)
                  PVIR(MA,MC) = PVIR(MA,MC) +
     *                          TERM*XPQRJ(MK+JSHIFT,MAI,MB+NOB)/DIKAB
                  WVIR(MA,MC) = WVIR(MA,MC) -
     *                          TERM*XPQRJ(MK+JSHIFT,MAI,MB+NOB)
                  XLAI1(MA,MI+NACORE)= XLAI1(MA,MI+NACORE) +
     *                          TERM*XPQRJ(MK+JSHIFT,MAC,MB+NOB)
  330          CONTINUE
  334       CONTINUE
  338    CONTINUE
  340 CONTINUE
C
  350 CONTINUE
C
C SYMMETRIZE THE CORE ORBITALS
C
      IF(NACORE.GT.0) THEN
         DO 390 MI=1,NOA
            DO 400 MJ=1,NACORE
                POCC(MI,MJ) = POCC(MJ,MI)
  400       CONTINUE
  390    CONTINUE
      END IF
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UMPGEN
      SUBROUTINE UMPGEN(E2,EA,EB,POCCA,POCCB,WOCCA,WOCCB,PVIRA,PVIRB,
     *                  WVIRA,WVIRB,NOA,NOB,NVIRA,NVIRB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00)
C
      DIMENSION EA(*),EB(*),POCCA(NOA,NOA),POCCB(NOB,NOB),
     *          WOCCA(NOA,NOA),WOCCB(NOB,NOB),PVIRA(NVIRA,NVIRA),
     *          PVIRB(NVIRB,NVIRB),WVIRA(NVIRA,NVIRA),WVIRB(NVIRB,NVIRB)
C
C     ----- UMP2 ENERGY (ROUTINE MODELED AFTER MP2GEN) -----
C
      E2 = ZERO
      DO 10 MI=1,NOA
         E2 = E2 + WOCCA(MI,MI)
   10 CONTINUE
      DO 15 MI=1,NOB
         E2 = E2 + WOCCB(MI,MI)
   15 CONTINUE
      E2 = -E2*PT5
C
C     ----- ENERGY CONTRIBUTIONS TO WOCCA,WOCCB,WVIRA,AND WVIRB -----
C
      DO 22 MJ=1,NOA
         DO 20 MI=1,NOA
            WOCCA(MI,MJ) = WOCCA(MI,MJ)
     *                     - POCCA(MI,MJ)*(EA(MI)+EA(MJ))*PT5
   20    CONTINUE
   22 CONTINUE
C
      DO 28 MJ=1,NOB
         DO 26 MI=1,NOB
            WOCCB(MI,MJ) = WOCCB(MI,MJ)
     *                     - POCCB(MI,MJ)*(EB(MI)+EB(MJ))*PT5
   26    CONTINUE
   28 CONTINUE
C
      DO 32 MB=1,NVIRA
         DO 30 MA=1,NVIRA
            WVIRA(MA,MB) = WVIRA(MA,MB) -
     *                     PVIRA(MA,MB)*(EA(MA+NOA)+EA(MB+NOA))*PT5
   30    CONTINUE
   32 CONTINUE
C
      DO 38 MB=1,NVIRB
         DO 36 MA=1,NVIRB
            WVIRB(MA,MB) = WVIRB(MA,MB) -
     *                     PVIRB(MA,MB)*(EB(MA+NOB)+EB(MB+NOB))*PT5
   36    CONTINUE
   38 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UMPGRD
      SUBROUTINE UMPGRD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000, MXRT=100, MXAO=8192)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,GOPARR,DSKWRK,MASWRK,
     *        ABEL,ABELPT,DIRSCF,FDIFF,DIRTRF
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      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 /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
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWORD,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NO,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFMEM/ IFCK,IDUM,IDEN,JDUM(3),IBUF,KDUM(3)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
C     COSMO INFORMATION
C
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
      LOGICAL ISEPS, USEPS
      DIMENSION QVPOT(NPPA)
      COMMON /ISEPS/  ISEPS, USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     1                DEN(MAXDEN),AR(LENABC),
     3                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /SOLVI / IATSP(LENABC+1),N0(2),NP1,NP2,ISKIP
      COMMON /COSMO1/ SE2,SECORR,ETOTS,CDUM,QVCOSMO,
     *                CSPOT(NPPA),ICORR,ITRIPO,ITRIP2,ITRIP3,ITRIP4,
     *                NATCOS,NQS,ITERC
      COMMON /COSMO2/ QENUC,ELAST,EMP2COS,EMP2COS2,ETOTSMP,SAVESE,
     *                EMP2LAST,MP2TRIP,MP2ITER,MP2FACT
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      DATA KRECA1,KRECA2,KRECA3,KRECA4,KRECA5,KRECA6
     *        /411,412,413,414,415,416/
      DATA KRECB1,KRECB2,KRECB3,KRECB4,KRECB5,KRECB6
     *        /421,422,423,424,425,426/
      DATA KRECA7,KRECA8,KRECA9/417,418,419/
      DATA KRECB7,KRECB8,KRECB9/427,428,429/
      DATA NFTAA,NFTAB,NFTBA,NFTBB,NFTIJ,NFT16/9,41,42,43,15,16/
C
      DATA CHECK/8HCHECK   /
C
C     ----- THIS CODE IS WRITTEN BY CHRISTINE AIKENS -----
C            IT IS MODELED AFTER THE MP2GRD ROUTINE.
C
C     ----- THIS ROUTINE COMPLETES THE SETUP FOR THE UMP2 GRADIENT -----
C     SOME OF THE 1E- DENSITY MATRICES AND THE ENERGY WERE CALCULATED
C     ALREADY (SEE UMPGE2)
C     FIRST, FINISH THE MP2 LAGRANGIAN AND SOLVE THE CPHF-LIKE EQUATION.
C     THEN, FINISH THE DENSITY MATRIX AND THE ENERGY-WEIGHTED DENSITY
C     MATRIX. AT THIS POINT, IT IS POSSIBLE TO DO 1E- PROPERTIES.
C     FINALLY, CREATE THE NON-SEPARABLE PART OF THE 2ND ORDER DENSITY.
C     THE SEPARABLE PART WILL BE DONE IN DABUMP, WHICH IS CALLED FROM
C     HFGRAD.
C
      SOME = MASWRK  .AND.  NPRINT.NE.-5
      IF(SOME) WRITE(IW,9992)
C
C        ---- OPTION TO DEFEAT SYMMETRY USAGE ----
C
      ABEL=ABELPT()
      IF(NOSYMT.EQ.1) CALL SYMOFF
C
      CUTOFF=TOL
C
      NOCA  = NOA - NACORE
      NOCB  = NOB - NBCORE
      NVIRA = NO  - NOA
      NVIRB = NO  - NOB
C
      NBF2  = (NBF*NBF+NBF)/2
      NBF3  = NBF*NBF
      NROTA = NVIRA*NOA
      NROTB = NVIRB*NOB
C
C     ----- OBTAIN ORBITAL SYMMETRY LABELS -----
C
C     IVECA:   NBF*NBF              - ALPHA ORBITAL COEFF
C     IVECB:   NBF*NBF              - BETA ORBITAL COEFF
C     IIRPA:   NBF                  - IRREP LABELS FOR ALPHA MO'S
C     IIRPB:   NBF                  - IRREP LABELS FOR BETA MO'S
C     IMOLAB:  NBF                  - TEMPORARY DATA
C     IMODEG:  NBF                  - TEMPORARY DATA
C     IS:      NBF2                 - OVERLAP INTEGRALS
C     IQ:      NBF2                 - Q MATRIX
C     IT:      NBF                  - TEMPORARY DATA
C
      CALL VALFM(LOADFM)
      IVECA = 1     + LOADFM
      IVECB = IVECA + NBF3
      IIRPA = IVECB + NBF3
      IIRPB = IIRPA + NBF
      LAST  = IIRPB + NBF
      NEEDA = LAST  - LOADFM
      CALL GETFM(NEEDA)
C
      CALL VALFM(LOADFM)
      IMOLAB = 1      + LOADFM
      IMODEG = IMOLAB + NBF
      IS     = IMODEG + NBF
      IQ     = IS     + NBF2
      IT     = IQ     + NBF3
      LAST   = IT     + NBF
      NEEDB  = LAST   - LOADFM
      CALL GETFM(NEEDB)
      IF(EXETYP.EQ.CHECK) GO TO 90
C
      CALL DAREAD(IDAF,IODA,X(IVECA),NBF3,15,0)
      CALL DAREAD(IDAF,IODA,X(IS   ),NBF2,12,0)
      CALL DAREAD(IDAF,IODA,X(IQ   ),NBF3,45,0)
      CALL TRFSYM(X(IMOLAB),X(IIRPA),X(IMODEG),X(IQ),X(IS),X(IVECA),
     *            X(IT),IA,NQMT,NBF,NO,NBF)
      CALL DAREAD(IDAF,IODA,X(IVECB),NBF3,19,0)
      CALL DAREAD(IDAF,IODA,X(IQ   ),NBF3,45,0)
      CALL TRFSYM(X(IMOLAB),X(IIRPB),X(IMODEG),X(IQ),X(IS),X(IVECB),
     *            X(IT),IA,NQMT,NBF,NO,NBF)
C
   90 CONTINUE
      CALL RETFM(NEEDB)
C
C
      NTSAVE = NT
      IF(DIRSCF  .AND.  .NOT.ABEL) NT=1
C
C     ----- NOW COMPLETE MP2 LAGRANGIAN -----
C
C     ALREADY DEFINED:
C     -IVECA-,-IVECB-   NBF*NBF              - ORBITAL COEFF     -
C     -IIRPA-,-IIRPB-   NBF                  - MOSIRP            -
C
C     NEED TO DEFINE:
C     -ILAIA-,-ILAIB-   NROTA,NROTB          - XLAI ALPHA,BETA   -
C     -IENGA-,-IENGB-   NBF                  - ORBITAL ENERGIES  -
C
C     ALSO NEED (TEMPORARILY):
C     -IRHOA-  MAX(NOA*NOA,NVIRA*NVIRA,NBF2) - PIJ,PAB,PMN,XLAI2 -
C     -IRHOB-  MAX(NOB*NOB,NVIRB*NVIRB,NBF2) - PIJ,PAB,PMN,XLAI2 -
C     -IXLMNA- NBF3                          - X,XLMNA           -
C     -IXLMNB- NBF3                          -   XLMNB           -
C     -IXX-    MAX(NBF3,NINTMX)              - XX,YA             -
C     -IYB-    NBF3                          -    YB             -
C     -IIX-    NINTMX                        - IX                -
C
      CALL VALFM(LOADFM)
      ILAIA = 1     + LOADFM
      ILAIB = ILAIA + MAX(NINTMX,NROTA)
      IENGA = ILAIB + MAX(NINTMX,NROTB)
      IENGB = IENGA + NBF
      LAST  = IENGB + NBF
      NEEDC = LAST  - LOADFM
      CALL GETFM(NEEDC)
C
      CALL VALFM(LOADFM)
      IRHOA = 1     + LOADFM
      IRHOB = IRHOA + MAX(NOA*NOA,NVIRA*NVIRA,NBF2)
      IXLMNA= IRHOB + MAX(NOB*NOB,NVIRB*NVIRB,NBF2)
      IXLMNB= IXLMNA+ NBF3
      IXX   = IXLMNB+ NBF3
      IIX   = IXX   + MAX(NBF3,NINTMX)
      IYB   = IIX   + NINTMX
      LAST  = IYB   + NBF3
      NEEDD = LAST  - LOADFM
      CALL GETFM(NEEDD)
      CALL TSECND(TIM0)
      IF(EXETYP.EQ.CHECK) GO TO 250
C
C     READ IN ORBITAL ENERGIES
C
      CALL DAREAD(IDAF,IODA,X(IENGA),NBF,17,0)
      CALL DAREAD(IDAF,IODA,X(IENGB),NBF,21,0)
C
      CALL UMPLAI(X(IRHOA),X(IRHOA),X(IRHOA),X(IRHOB),X(IRHOB),X(IRHOB),
     *            X(IXLMNA),X(IXX),X(IYB),X(IXLMNA),X(IXLMNB),X(IXX),
     *            X(IIX),X(ILAIA),X(IRHOA),X(ILAIB),X(IRHOB),
     *            NOA,NVIRA,NOB,NVIRB,NBF,KRECA1,KRECA2,KRECA3,KRECA4,
     *            KRECB1,KRECB2,KRECB3,KRECB4,X(IVECA),X(IVECB))
C
      IF(DBUG) THEN
         IF(SOME) WRITE(IW,9999)
         CALL PRSQ(X(ILAIA),NOA,NVIRA,NVIRA)
         IF(SOME) WRITE(IW,9899)
         CALL PRSQ(X(ILAIB),NOB,NVIRB,NVIRB)
      END IF
C
  250 CONTINUE
      CALL TSECND(TIM1)
      TIME = TIM1 - TIM0
      IF(SOME) WRITE(IW,9997) TIME
      CALL FLSHBF(IW)
C
      CALL RETFM(NEEDD)
C
C     ----- MP2 CPHF EQUATIONS (AO BASIS) -----
C
      WRITE(IW,9980)
      CALL VALFM(LOADFM)
      IRHOA = 1     + LOADFM
      IRHOB = IRHOA + NBF2
      IXX   = IRHOB + NBF2
      IYA   = IXX   + NBF3
      IYB   = IYA   + NBF3
      IZA   = IYB   + NBF3
      IZB   = IZA   + NBF3
      IXLMNA= IZB   + NBF3
      IXLMNB= IXLMNA+ NBF3
      LBUF  = IXLMNB+ NBF3
      IBUF  = LBUF  + NINTMX
      IPAI1 = IBUF  + NINTMX
      IPAI2 = IPAI1 + NROTA
      IPBJ1 = IPAI2 + NROTB
      IPBJ2 = IPBJ1 + NROTA
      IERR  = IPBJ2 + NROTB
      LAST  = IERR  + NBF3
      NEEDE = LAST  - LOADFM
      CALL GETFM(NEEDE)
C
      MAXIO  = 3*MAXDII
      ITERMX = MAXITC
      ITERMX2= (ITERMX*ITERMX+ITERMX)/2
C
      CALL VALFM(LOADFM)
      IWRK   = 1     + LOADFM
      IAM    = IWRK  + NBF3
      IX     = IAM   + MAXDII*MAXDII
      IPVT   = IX    + ITERMX
      IB     = IPVT  + ITERMX
      IODIIS = IB    + ITERMX2
      LAST   = IODIIS+ MAXIO
      NEEDF  = LAST  - LOADFM
      CALL GETFM(NEEDF)
C
      IF(EXETYP.EQ.CHECK) GO TO 350
C
      CALL DAREAD(IDAF,IODA,X(IENGA),NBF,17,0)
      CALL DAREAD(IDAF,IODA,X(IENGB),NBF,21,0)
C
      CALL CPHFU(X(IRHOA),X(IRHOB),X(IYA),X(IYB),X(IXLMNA),X(IXLMNB),
     *           X(IXX),X(LBUF),X(IBUF),X(ILAIA),X(IPAI1),X(IPBJ1),
     *           X(ILAIB),X(IPAI2),X(IPBJ2),NOA,NVIRA,NOB,NVIRB,NROTA,
     *           NROTB,X(IVECA),X(IENGA),X(IVECB),X(IENGB),X(IWRK),
     *           X(IAM),X(IX),X(IB),X(IPVT),
     *           X(IODIIS),NBF,X(IERR),ITERMX,X(IZA),X(IZB))
C
      CALL DCOPY(NROTA,X(IPAI1),1,X(ILAIA),1)
      CALL DCOPY(NROTB,X(IPAI2),1,X(ILAIB),1)
C
  350 CONTINUE
      CALL RETFM(NEEDF)
      CALL RETFM(NEEDE)
C
      NT = NTSAVE
      IF(NOSYMT.EQ.1) CALL SYMON
C
      CALL TSECND(TIM2)
      TIME = TIM2 - TIM1
      IF(SOME) WRITE(IW,9996) TIME
      CALL FLSHBF(IW)
C
C     ----- COMPLETE W2 -----
C
C     FOR ALPHA:
C     ILAIA =  NROTA        PAI, OBTAINED FROM CPHF EQUATION
C     IENGA =  NBF          ORBITAL ENERGY
C     IPOCCA=  NOA*NOA      OCCUPIED-OCCUPIED DENSITY (POCCA)
C     IPVIRA=  NVIRA*NVIRA  VIRTUAL-VIRTUAL DENSITY (PVIRA)
C     IWOCCA=  NOA*NOA      ENERGY WEIGHT DENSITY (WOCCA)
C     IWAIA =  NOA*NVIRA    ENERGY WEIGHT DENSITY (WAIA)
C     PLUS CORRESPONDING BETA AND:
C     IXX   =  NBF2         XX
C     IX    =  NBF3         X
C     IWRK1 =  NBF3         SCRATCH
C     IWRK2 =  NBF3         SCRATCH
C
      CALL VALFM(LOADFM)
      IPOCCA  = 1     + LOADFM
      IPOCCB  = IPOCCA+ NOA*NOA
      IPVIRA  = IPOCCB+ NOB*NOB
      IPVIRB  = IPVIRA+ NVIRA*NVIRA
      IWOCCA  = IPVIRB+ NVIRB*NVIRB
      IWOCCB  = IWOCCA+ NOA*NOA
      IWAIA   = IWOCCB+ NOB*NOB
      IWAIB   = IWAIA + NVIRA*NOA
      IRHOA   = IWAIB + NVIRB*NOB
      IRHOB   = IRHOA + NBF2
      IXA     = IRHOB + NBF2
      IXB     = IXA   + NBF3
      IWRK1   = IXB   + NBF3
      IWRK2   = IWRK1 + NBF3
      IWRK3   = IWRK2 + NBF3
      NBUF    = IWRK3 + NBF3
      IBUF    = NBUF + NINTMX
      LAST    = IBUF + NINTMX
      NEEDG   = LAST  - LOADFM
      CALL GETFM(NEEDG)
      IF(EXETYP.EQ.CHECK) GO TO 550
C
      CALL DAREAD(IDAF,IODA,X(IPOCCA),NOA  *NOA  ,KRECA1,0)
      CALL DAREAD(IDAF,IODA,X(IPVIRA),NVIRA*NVIRA,KRECA2,0)
      CALL DAREAD(IDAF,IODA,X(IWAIA ),NVIRA*NOA  ,KRECA3,0)
      CALL DAREAD(IDAF,IODA,X(IWOCCA),NOA  *NOA  ,KRECA5,0)
      CALL DAREAD(IDAF,IODA,X(IPOCCB),NOB  *NOB  ,KRECB1,0)
      CALL DAREAD(IDAF,IODA,X(IPVIRB),NVIRB*NVIRB,KRECB2,0)
      CALL DAREAD(IDAF,IODA,X(IWAIB ),NVIRB*NOB  ,KRECB3,0)
      CALL DAREAD(IDAF,IODA,X(IWOCCB),NOB  *NOB  ,KRECB5,0)
C
C     ----- TRANSFORM PIJ, PAB, AND PAI TO P2(MUNU) -----
C
      CALL UMPGTR(X(IPOCCA),X(IPVIRA),X(ILAIA),X(IVECA),NBF,NOA,NVIRA,
     *            X(IWRK1),X(IWRK2),X(IRHOA))
      CALL DAWRIT(IDAF,IODA,X(IRHOA),NBF2,KRECA7,0)
      CALL UMPGTR(X(IPOCCB),X(IPVIRB),X(ILAIB),X(IVECB),NBF,NOB,NVIRB,
     *            X(IWRK1),X(IWRK2),X(IRHOB))
      CALL DAWRIT(IDAF,IODA,X(IRHOB),NBF2,KRECB7,0)
C
C     DO PART 3 OF WOCC AND PART 2 OF WAI
C
      CALL FNSHW2(X(IPOCCA),X(IPOCCB),X(ILAIA),X(ILAIB),X(IWOCCA),
     *            X(IWOCCB),X(IWAIA),X(IWAIB),X(IENGA),X(IENGB),
     *            X(IRHOA),X(IRHOB),X(IXA),X(IXB),X(IWRK1),
     *            X(IWRK2),X(IWRK3),X(NBUF),X(IBUF),X(IVECA),X(IVECB),
     *            NOA,NOB,NVIRA,NVIRB,NBF,NBF2)
C
C     NEXT, READ IN WVIR
C
      CALL DAREAD(IDAF,IODA,X(IXA),NVIRA*NVIRA,KRECA6,0)
      CALL DAREAD(IDAF,IODA,X(IXB),NVIRB*NVIRB,KRECB6,0)
C
C     ----- TRANSFORM WIJ, WAB, AND WAI TO P2(MUNU) -----
C
      CALL UMPGTR(X(IWOCCA),X(IXA),X(IWAIA),X(IVECA),NBF,NOA,NVIRA,
     *            X(IWRK1),X(IWRK2),X(IRHOA))
      CALL UMPGTR(X(IWOCCB),X(IXB),X(IWAIB),X(IVECB),NBF,NOB,NVIRB,
     *            X(IWRK1),X(IWRK2),X(IRHOB))
      CALL UMPW2(X(IVECA),X(IENGA),X(IWRK1),X(IRHOA),NOA,NBF,
     *           KRECA7,KRECA8,KRECA9,16)
      CALL UMPW2(X(IVECB),X(IENGB),X(IWRK1),X(IRHOB),NOB,NBF,
     *           KRECB7,KRECB8,KRECB9,20)
C
C         SPIN CONTAMINATION CALCULATION (USES OVERLAP INTEGRALS)
C         VALUES FROM RELAXED DENSITY ARE HIGHER THAN THE HF VALUES.
C         IT WOULD BE MORE BETTER TO USE EXPECTATION VALUE DENSITIES.
C
C     CALL DAREAD(IDAF,IODA,X(IWRK1),NBF2,12,0)
C     CALL DAREAD(IDAF,IODA,X(IRHOA),NBF2,16,0)
C     CALL DAREAD(IDAF,IODA,X(IRHOB),NBF2,20,0)
C     CALL SPIN(SZ,S2,X(IRHOA),X(IRHOB),X(IWRK1),X(IWRK2),X(IWRK3),
C    *          IA,NOA,NOB,NBF,NBF2)
C     IF(MASWRK) WRITE(IW,9910) S2
C9910 FORMAT(1X,'THE FIRST ORDER WAVEFUNCTION HAS <S**2>=',F8.3)
C
      CALL TSECND(TIM3)
      TIME = TIM3 - TIM2
      IF(SOME) WRITE(IW,9995) TIME
      IF(MASWRK)
     *    WRITE(IW,*) '..... DONE WITH UMP2 CPHF AND DENSITY .....'
      CALL TIMIT(1)
      CALL FLSHBF(IW)
C
  550 CONTINUE
      CALL RETFM(NEEDG)
C
      IF(NPRINT.NE.-5) CALL MP2NOS
C
      ETOT = EMP2
      CALL DAWRIT(IDAF,IODA,ENUC,MXRT+15,2,0)
C
C     COSMO CHANGES - KIM, PUT IN BY VOLKER NOV.1998
C
C     COMPUTE PROPERTIES INCLUDING THE COSMO SOLVATION CORRECTION
C     FOR THE SECOND ORDER WAVE FUNCTION AND ITERATE BOTH Q'S
C     AND UMP2 ENERGY TO SELF-CONSISTENCY
C     NOTE: THIS PART OF THE ROUTINE IS DIRECTLY FROM MP2GRD
C
      IF(ISEPS) THEN
          WRITE(IW,*) "MP2ITER IN MP2GRD = ",MP2ITER
          WRITE(IW,*)
          WRITE(IW,*)"AMOS'S DISTRIBUTED MULTIPOLE ANALYSIS:"
          WRITE(IW,*)"      -- CALLING FROM UMPGRD --"
          WRITE(IW,*)
C
          CALL NUMPROP
C
          WRITE(IW,*)"NEW SURFACE CHARGES/SCREENING E FOR UMP2",MP2TRIP
          IF(MP2TRIP.EQ.1) THEN
             NDEN2=(NDEN*NDEN+NDEN)/2
             LENAB2 = NPS2 + NDEN*NPS + NDEN2
             CALL VALFM(LOADFM)
             LABCMT = LOADFM + 1
             LAST   = LABCMT + LENAB2
             NEEDCO = LAST - LOADFM - 1
             CALL GETFM(NEEDCO)
C
C              RECOVER PREVIOUS A,B,C MATRICES
C
             NFT26 = 26
             LENB = NDEN*NPS
             CALL SQREAD(NFT26,X(LABCMT)              ,NPS2)
             CALL SQREAD(NFT26,X(LABCMT+NPS2)         ,LENB)
             CALL SQREAD(NFT26,X(LABCMT+NPS2+NDEN*NPS),NDEN2)
             CALL SEQREW(NFT26)
C
C              AND COMPUTE THE SOLVATION ENERGY WITH THEM
C
             CALL DECORR(X(LABCMT),LENAB2)
             CALL RETFM(NEEDCO)
C
             DSECORR=SECORR-SAVESE
C            EMP2SOL1=(EMP2-ETOTSMP)*627.517D+00
C            EMP2SOL=(EMP2-ETOTSMP-SECORR)*627.517D+00
C            EMP2SOL2=(EMP2-ETOTSMP-DSECORR)*627.517D+00
      CALL FLSHBF(6)
C
C   NOW WE HAVE TO DO THE OCE CORRECTION ANALOGOUS TO THE SCF PART
C   QVCOSMO SHOULD BE THERE FROM DECORR, WE NEED QVGAMESS FROM THE
C   MP2 DENSITY - THIS IS FOR THE CONVERGED MP2 RESULTS, MP2ITER=0
C
             EOC1=0.0D+00
             IF(MP2ITER.EQ.0) THEN
                DO 5 II=1,NPS
                   QVPOT(II)=0.0D+00
  5             CONTINUE
      CALL FLSHBF(6)
C     CALL VALFM(LOADFM)
C     LAST  = 1   + LOADFM
C     LAST  = LAST + NQS
C     NEEDG2= LAST - LOADFM
C     CALL GETFM(NEEDG2)
                CALL COSPOT(QVPOT)
C
C   COSPOT GAVE US V'; NOW FORM QV'
C
                QVCORR=0.0D+00
                DO 6 I=1,NQS
                   QVCORR = QVCORR + QVPOT(I)*COSZAN(I)
  6             CONTINUE
                WRITE(IW,*)"QV (GAMESS,UMP2): QVCORR = ",QVCORR
                WRITE(IW,*)"QV (COSMO,UMP2) : QVCOSMO= ",QVCOSMO
                EOC1=-QVCORR+QVCOSMO
                EOC2=EOC1*627.517D+00
                CHECK2=QVCORR*QVCORR-QVCOSMO*QVCOSMO
                WRITE(IW,*)"-QV(GAMESS,UMP2)+QV(COSMO,UMP2):",EOC1
                WRITE(IW,*)"POTENTIAL NORM DIFFERENCE :",CHECK2
                WRITE(IW,*)"QVCOSMO/QVGAMESS:",QVCOSMO/QVCORR
                WRITE(IW,*)
                WRITE(IW,9400)
                WRITE(IW,*)"  OUTLYING CHARGE EFFECT (A.U.)      =",EOC1
                WRITE(IW,*)"  OUTLYING CHARGE EFFECT (KCAL/MOL)  =",EOC2
                WRITE(IW,9400)
             ENDIF
      CALL FLSHBF(6)
C
C   EOC1 IS NOW THE MP2 OUTLYING CHARGE CORRECTION, WHICH MUST BE
C   ADDED TO EMP2 AFTER PRINTOUT
C
             EMP2C=EMP2+EOC1
             WRITE(IW,9500)ESCF,SAVESE,SECORR,DSECORR,EMP2,EMP2C
             EMP2=EMP2C
         ENDIF
C
         IF(MP2ITER.EQ.1) THEN
             WRITE(IW,*)" 1 ELECTRON CONTRIBUTION FROM MP2 CHARGES:"
             CALL ONEEI
         ENDIF
C
         WRITE(IW,*)
         WRITE(IW,*) "**********COSMO INFORMATION**UMPGRD**************"
         WRITE(IW,*)
         WRITE(IW,*)"  NUMBER OF SURFACE SEGMENTS IS ",NPS
         WRITE(IW,*)
         AREA=0.0D+00
         DO 889 I=1,NPS
           AREA=AREA+AR(I)
           WRITE(IW,'(2I4,6F10.5)') I,IATSP(I),CORZAN(1,I),CORZAN(2,I),
     *CORZAN(3,I),COSZAN(I),AR(I),COSZAN(I)/AR(I)
 889       CONTINUE
           WRITE(IW,*)
           WRITE(IW,*)" TOTAL SURFACE AREA OF CAVITY (A**2)     = ",AREA
           WRITE(IW,*)
           WRITE(IW,*)"***********************************************"
           WRITE(IW,*)
C
C     MP2ITER=1 MEANS THAT THE UMP2 ENERGY IS NOT CONVERGED, SO WE GO
C     BACK HERE AFTER THE UMP2 DENSITY IS UPDATED
C
         WRITE(IW,*)"RETURNING FROM UMPGRD FOR MP2ITER=1",MP2ITER
         IF(MP2ITER.EQ.1) RETURN
      END IF   ! END OF COSMO INSERTION
C
C        A RUN TO GENERATE ONLY THE DENSITY CAN NOW EXIT.
C
      IF(MPPROP.EQ.1) THEN
         CALL RETFM(NEEDC)
         CALL RETFM(NEEDA)
         RETURN
      END IF
C
C     ----- NON-SEPARABLE CONTRIBUTIONS TO THE 2ND ORDER DENSITY -----
C     IAIBJ  = NVIR*NOC*NVIR       AIBJ
C     IWRK1  = NBF3                X
C     IWRK2  = NBF3                Y,XX
C     IXIN   = NBF2*NBF            XIN
C
      CALL VALFM(LOADFM)
      IAIBJ = 1   + LOADFM
      IAIBJ2= IAIBJ + MAX(NVIRA*NOCA*NVIRA,NVIRB*NOCB*NVIRB)
      IWRK1 = IAIBJ2+ MAX(NVIRB*NOCA*NVIRA,NVIRA*NOCB*NVIRB)
      IWRK2 = IWRK1 + NBF3
      IXIN  = IWRK2 + NBF3
      IENGA = IXIN  + NBF2*NBF
      IENGB = IENGA + NBF
      LAST  = IENGB + NBF
      NEEDG = LAST - LOADFM
C
      CALL GOTFM(NGOTMX)
      IF(NEEDG.GT.NGOTMX) THEN
         IF(MASWRK) WRITE(IW,9993) NEEDG-NGOTMX
         CALL ABRT
      END IF
      CALL GETFM(NEEDG)
      IF(EXETYP.EQ.CHECK) GO TO 750
C
      CALL DAREAD(IDAF,IODA,X(IENGA),NBF,17,0)
      CALL DAREAD(IDAF,IODA,X(IENGB),NBF,21,0)
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTIJ)
C
C     FORM PARTIALLY TRANSFORMED NON-SEPARABLE CONTRIBUTION
C
      IF (SOME) WRITE(IW,9990)
      CALL NONSEP(X(IAIBJ),X(IAIBJ2),X(IENGA),X(IENGB),X(IVECA),
     *           X(IVECB),X(IWRK1),
     *           X(IWRK2),X(IWRK2),X(IXIN),NBF,NBF2,NOA,NOB,NVIRA,NVIRB,
     *           NFTAA,NFTAB,NFT16,CUTOFF,SOME,NACORE,NBCORE,NOCA,NOCB)
      CALL NONSEP(X(IAIBJ),X(IAIBJ2),X(IENGB),X(IENGA),X(IVECB),
     *           X(IVECA),X(IWRK1),
     *           X(IWRK2),X(IWRK2),X(IXIN),NBF,NBF2,NOB,NOA,NVIRB,NVIRA,
     *           NFTBB,NFTBA,NFTIJ,CUTOFF,SOME,NBCORE,NACORE,NOCB,NOCA)
C
  750 CONTINUE
      CALL RETFM(NEEDG)
      CALL RETFM(NEEDC)
      CALL RETFM(NEEDA)
C
      IF(MASWRK) WRITE(IW,*) '..... DONE WITH MP2 2DM .....'
      CALL TIMIT(1)
      RETURN
C
 9999 FORMAT(/' --- ALPHA LAGRANGIAN ---'/)
 9899 FORMAT(/' --- BETA  LAGRANGIAN ---'/)
 9997 FORMAT(/' -UMP2- LAGRANGIAN     TOOK = ',F10.3,' SECONDS')
 9996 FORMAT( ' -UMP2- CPHF EQUATIONS TOOK = ',F10.3,' SECONDS')
 9995 FORMAT( ' -UMP2- W2(OCC,OCC)    TOOK = ',F10.3,' SECONDS')
 9993 FORMAT(/' OUT OF MEMORY FOR 2DM TRANSFORMATION'/
     *        ' INCREASE MEMORY SIZE BY AT LEAST ',I10,' WORDS'/)
 9992 FORMAT(/10X,21(1H-)/10X,'CPHF FOR UMP2 GRADIENT',
     *       /10X,21(1H-))
 9990 FORMAT(/1X,'FORMING NON-SEPARABLE PART OF THE UMP2 2ND ORDER',
     *           ' DENSITY MATRIX...')
 9980 FORMAT(/' ..... SOLVING CPHF IN THE AO BASIS .....')
 9400 FORMAT(2X,'***************************************************',
     *'**************')
 9500 FORMAT(/1X,'RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE'/
     *       8X,'  ECOSMO(0)                      =',1X,F20.10/
     *       8X,'  UHF SECORR                     =',1X,F20.10/
     *       8X,'  MP2 SECORR                     =',1X,F20.10/
     *       8X,'DELTA SECORR                     =',1X,F20.10/
     *       8X,'E(COSMO-MP2)                     =',1X,F20.10/
     *       8X,'E(COSMO-MP2,OCE-CORRECTED)       =',1X,F20.10/)
      END
C
C*MODULE MP2GRD  *DECK UMPLAI
      SUBROUTINE UMPLAI(PIJ1,PAB1,PMNA,PIJ2,PAB2,PMNB,X,YA,YB,XLMNA,
     *                  XLMNB,BUF,IBUF,XLAI1A,XLAI2A,XLAI1B,XLAI2B,
     *                  NOA,NVIRA,NOB,NVIRB,NBF,
     *                  KRECA1,KRECA2,KRECA3,KRECA4,
     *                  KRECB1,KRECB2,KRECB3,KRECB4,CA,CB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E
C
      DIMENSION PIJ1(NOA,NOA),PAB1(NVIRA,NVIRA),PMNA(*),CA(NBF,*)
      DIMENSION PIJ2(NOB,NOB),PAB2(NVIRB,NVIRB),PMNB(*),CB(NBF,*)
      DIMENSION X(NBF,NBF),XLMNA(*),XLMNB(*)
      DIMENSION YA(NBF,NBF),YB(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
      DIMENSION XLAI1A(NVIRA,NOA),XLAI2A(NVIRA,NOA)
      DIMENSION XLAI1B(NVIRB,NOB),XLAI2B(NVIRB,NOB)
C
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ONE=1.0D+00,HALF=0.5D+00)
C
C     ----- COMPLETE CONSTRUCTION OF THE MP2 LAGRANGIAN -----
C     ACCORDING TO EQUATION 11 OF THE FRISCH PAPER.
C
C        BEGIN TRANSFORMING -PIJ- AND -PAB- TO AO BASIS AS -PMN-
C
C     DO ALPHA FIRST:
C
      CALL DAREAD(IDAF,IODA,PIJ1,NOA*NOA,KRECA1,0)
      CALL MRARTR(PIJ1,NOA,NOA,NOA,CA,NBF,NBF,X,NBF)
      CALL MRARBR(CA,NBF,NBF,NOA,X,NBF,NBF,YA,NBF)
C
      CALL DAREAD(IDAF,IODA,PAB1,NVIRA*NVIRA,KRECA2,0)
      CALL MRARTR(PAB1,NVIRA,NVIRA,NVIRA,CA(1,NOA+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRA,ONE,CA(1,NOA+1),NBF,X,NBF,
     *           ONE,YA,NBF)
C
C     NOW BETA:
C
      CALL DAREAD(IDAF,IODA,PIJ2,NOB*NOB,KRECB1,0)
      CALL MRARTR(PIJ2,NOB,NOB,NOB,CB,NBF,NBF,X,NBF)
      CALL MRARBR(CB,NBF,NBF,NOB,X,NBF,NBF,YB,NBF)
C
      CALL DAREAD(IDAF,IODA,PAB2,NVIRB*NVIRB,KRECB2,0)
      CALL MRARTR(PAB2,NVIRB,NVIRB,NVIRB,CB(1,NOB+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRB,ONE,CB(1,NOB+1),NBF,X,NBF,
     *           ONE,YB,NBF)
C
      MUNU = 0
      DO 110 MU=1,NBF
         DO 100 NU=1,MU
            MUNU=MUNU+1
            PMNA(MUNU) = HALF*(YA(MU,NU) + YA(NU,MU))
            PMNB(MUNU) = HALF*(YB(MU,NU) + YB(NU,MU))
  100    CONTINUE
  110 CONTINUE
C
C     THIRD AND FOURTH TERM LAGRANGIAN CONTRIBUTIONS, IN AO BASIS
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL FLMATU(PMNA,YA,XLMNA,PMNB,YB,XLMNB,BUF,IBUF,NBF)
C
C     AFTER FLMATU, YA AND YB ARE THE SQUARE DENSITY MATRICES
C     FOR CLOSED SHELL, YA=YB=Y
C     ---- TRANSFORM THESE TWO TERMS TO MO BASIS,
C          AND ADD ALL CONTRIBUTIONS TOGETHER ----
C
C     FIRST ALPHA:
C
      CALL DAREAD(IDAF,IODA,XLAI1A,NVIRA*NOA,KRECA4,0)
C
      CALL MRARBR(YA,NBF,NBF,NBF,CA,NBF,NOA,X,NBF)
      CALL DGEMM('T','N',NVIRA,NOA,NBF,ONE,CA(1,NOA+1),NBF,X,NBF,
     *           ONE,XLAI1A,NVIRA)
C
      CALL DAREAD(IDAF,IODA,XLAI2A,NVIRA*NOA,KRECA3,0)
C
      DO 30 MI=1,NOA
      DO 30 MA=1,NVIRA
   30    XLAI1A(MA,MI) = -( XLAI1A(MA,MI)+XLAI2A(MA,MI) )
C
C     NOW BETA:
C
      CALL DAREAD(IDAF,IODA,XLAI1B,NVIRB*NOB,KRECB4,0)
C
      CALL MRARBR(YB,NBF,NBF,NBF,CB,NBF,NOB,X,NBF)
      CALL DGEMM('T','N',NVIRB,NOB,NBF,ONE,CB(1,NOB+1),NBF,X,NBF,
     *           ONE,XLAI1B,NVIRB)
C
      CALL DAREAD(IDAF,IODA,XLAI2B,NVIRB*NOB,KRECB3,0)
C
      DO 40 MI=1,NOB
      DO 40 MA=1,NVIRB
   40    XLAI1B(MA,MI) = -( XLAI1B(MA,MI)+XLAI2B(MA,MI) )
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK FLMATU
      SUBROUTINE FLMATU(PMNA,YA,XLMNA,PMNB,YB,XLMNB,BUF,IBUF,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E,SCHWRZ,DIRSCF,FDIFF
C
      DIMENSION PMNA(*),XLMNA(*),PMNB(*),XLMNB(*)
      DIMENSION YA(NBF,NBF),YB(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXAO=8192)
C
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      PARAMETER (HALF=0.5D+00)
C
      DATA UHF/8HUHF     /
C
C           NOW FORM FOCK-LIKE MATRIX
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      IF(DIRSCF) THEN
C
         SCHWRZ=ISCHWZ.GT.0
         NINT=0
         NSCHWZ=0
         L1 = NBF
         L2 = (NBF*NBF+NBF)/2
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
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
         CALL VALFM(LOADFM)
         LGHOND = LOADFM + 1
         LXINTS = LGHOND + MAXG
         LDSH   = LXINTS + NSH2
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
         IF(SCHWRZ) THEN
            DUMMY = 0.0D+00
            CALL SHLDEN(UHF,PMNA,PMNB,DUMMY,XX(LDSH),IA,L1,L2,NSH2,1)
            CALL DAREAD(IDAF,IODA,XX(LXINTS),NSH2,54,0)
         END IF
         CALL VCLR(XLMNA,1,L2)
         CALL VCLR(XLMNB,1,L2)
         CALL TWOEI(UHF,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     *              INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *              DUMMY,DUMMY,IDUMM,1,XX(LXINTS),NSH2,XX(LGHOND),
     *              MAXG,XX(LDDIJ),IA,PMNA,XLMNA,PMNB,XLMNB,XX(LDSH),
     *              DUMMY,DUMMY,1)
         CALL DSCAL(L2,HALF,XLMNA,1)
         CALL DSCAL(L2,HALF,XLMNB,1)
         II=0
         DO 210 I=1,L1
            II = II+I
            XLMNA(II) = XLMNA(II) + XLMNA(II)
            XLMNB(II) = XLMNB(II) + XLMNB(II)
  210    CONTINUE
C
         CALL RETFM(NEED)
      ELSE
        CALL SEQREW(IJK)
        CALL HSTARU(PMNA,XLMNA,PMNB,XLMNB,BUF,IBUF,YB,YB,NINTMX,IA,NOPK)
C
C   NOTE: YB IS JUST A FILLER.  WE DON'T REQUIRE XP AND XK SINCE NOPK=T
C
      END IF
C
C     ----- SYMMETRIZE FOCK-LIKE MATRIX -----
C
      CALL SYMH(XLMNA,YA,IA)
      CALL SYMH(XLMNB,YB,IA)
C
      CALL VCLR(YA,1,NBF*NBF)
      CALL VCLR(YB,1,NBF*NBF)
      MUNU = 0
      DO 20 MU=1,NBF
      DO 20 NU=1,MU
         MUNU=MUNU+1
         YA(MU,NU) = XLMNA(MUNU)
         YA(NU,MU) = XLMNA(MUNU)
         YB(MU,NU) = XLMNB(MUNU)
   20    YB(NU,MU) = XLMNB(MUNU)
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK CPHFU
      SUBROUTINE CPHFU(PMNA,PMNB,YA,YB,XLMNA,XLMNB,XX,BUF,IBUF,
     *                 XLAI1,PAI1,PBJ1,XLAI2,PAI2,PBJ2,NOA,NVIRA,NOB,
     *                 NVIRB,NROTA,NROTB,CA,EA,CB,EB,WRK,A,X,
     *                 B,IPVT,IODIIS,NBF,ERR,MAXIT,ZA,ZB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SOME,OUT,DBUG,PACK2E,NOTOPN
C
      DIMENSION PAI1(NVIRA,NOA),PBJ1(NVIRA,NOA),PMNA(*)
      DIMENSION PAI2(NVIRB,NOB),PBJ2(NVIRB,NOB),PMNB(*)
      DIMENSION XX(NBF,NBF),XLMNA(*),XLMNB(*)
      DIMENSION XLAI1(NVIRA,NOA),XLAI2(NVIRB,NOB)
      DIMENSION YA(NBF,NBF),YB(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
      DIMENSION ZA(NBF,NBF),ZB(NBF,NBF)
      DIMENSION CA(NBF,*),CB(NBF,*),EA(NBF),EB(NBF)
      DIMENSION ERR(NBF,NBF)
      DIMENSION WRK(*),A(*),X(*),IPVT(*),B(*)
      DIMENSION IODIIS(*)
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PRT/ SOME,OUT,DBUG
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, CVGTOL=1.0D-10)
      PARAMETER (HALF=0.5D+00)
C
      NBF3   = NBF*NBF
      ITDIIS = 1
      NOCP   = 0
      MAXIO  = 3*MAXDII
      MAXIT2 = (MAXIT*MAXIT+MAXIT)/2
      ERDIIS = ZERO
      NOTOPN = .TRUE.
C
C     SET UP AN INITIAL GUESS AT RESPONSE VECTORS (BOTH ALPHA AND BETA)
C
      CALL DCOPY(NROTA,XLAI1,1,PBJ1,1)
      DO 60 NI=1,NOA
         DO 50 NA=1,NVIRA
            PBJ1(NA,NI)=PBJ1(NA,NI)/(EA(NA+NOA)-EA(NI))
   50    CONTINUE
   60 CONTINUE
      CALL DCOPY(NROTB,XLAI2,1,PBJ2,1)
      DO 80 NI=1,NOB
         DO 70 NA=1,NVIRB
            PBJ2(NA,NI)=PBJ2(NA,NI)/(EB(NA+NOB)-EB(NI))
   70    CONTINUE
   80 CONTINUE
C
C     ----- FORM D(LAMDA,SIGMA) WHICH WE CALL PMN -----
C     ACCORDING TO EQUATION 17 OF THE FRISCH PAPER.
C
C        BEGIN TRANSFORMING -PBJ- TO AO BASIS AS -PMN-
C
      CALL VCLR(YA,1,NBF*NBF)
      CALL VCLR(YB,1,NBF*NBF)
C
      CALL MRARTR(PBJ1,NVIRA,NVIRA,NOA,CA,NBF,NBF,XX,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRA,TWO,CA(1,NOA+1),NBF,XX,NBF,
     *     ONE,YA,NBF)
      CALL MRARTR(PBJ2,NVIRB,NVIRB,NOB,CB,NBF,NBF,XX,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRB,TWO,CB(1,NOB+1),NBF,XX,NBF,
     *     ONE,YB,NBF)
C
C  START DIIS ITERATIONS
C
      DELTA = ZERO
      ITER   = 0
  100 CONTINUE
      ITER=ITER+1
      MUNU = 0
      DO MU=1,NBF
         DO NU=1,MU
            MUNU=MUNU+1
            PMNA(MUNU) = HALF*(YA(MU,NU) + YA(NU,MU))
            PMNB(MUNU) = HALF*(YB(MU,NU) + YB(NU,MU))
         ENDDO
      ENDDO
C
C           NOW FORM FOCK-LIKE MATRIX
C           DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C           STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL FLMATU(PMNA,ZA,XLMNA,PMNB,ZB,XLMNB,BUF,IBUF,NBF)
C
C     ---- TRANSFORM THE FOCK-LIKE MATRIX Z(MU,NU) TO MO BASIS ----
C
      CALL VCLR(PAI1,1,NROTA)
      CALL VCLR(PAI2,1,NROTB)
C
      CALL MRARBR(ZA,NBF,NBF,NBF,CA,NBF,NOA,XX,NBF)
      CALL DGEMM('T','N',NVIRA,NOA,NBF,ONE,CA(1,NOA+1),NBF,XX,NBF,
     *           ONE,PAI1,NVIRA)
C
      CALL MRARBR(ZB,NBF,NBF,NBF,CB,NBF,NOB,XX,NBF)
      CALL DGEMM('T','N',NVIRB,NOB,NBF,ONE,CB(1,NOB+1),NBF,XX,NBF,
     *           ONE,PAI2,NVIRB)
C
C     ---- CALCULATE PAI USING EQUATION 10 IN FRISCH'S PAPER ----
C
      DO 700 NI=1,NOA
         DO 710 NA=1,NVIRA
            PAI1(NA,NI)=(XLAI1(NA,NI)-PAI1(NA,NI))/(EA(NA+NOA)-EA(NI))
  710    CONTINUE
  700 CONTINUE
      DO 720 NI=1,NOB
         DO 730 NA=1,NVIRB
            PAI2(NA,NI)=(XLAI2(NA,NI)-PAI2(NA,NI))/(EB(NA+NOB)-EB(NI))
  730    CONTINUE
  720 CONTINUE
C
C     ---- TRANSFORM PAI ALPHA AND BETA TO AO BASIS ----
C
      CALL VCLR(ZA,1,NBF*NBF)
      CALL VCLR(ZB,1,NBF*NBF)
C
      CALL MRARTR(PAI1,NVIRA,NVIRA,NOA,CA,NBF,NBF,XX,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRA,TWO,CA(1,NOA+1),NBF,XX,NBF,
     *     ONE,ZA,NBF)
      CALL MRARTR(PAI2,NVIRB,NVIRB,NOB,CB,NBF,NBF,XX,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIRB,TWO,CB(1,NOB+1),NBF,XX,NBF,
     *     ONE,ZB,NBF)
C
      IF(ITER.GE.2) THEN
C
C     --- DETERMINE ALPHA ERROR VECTOR FOR DIIS ----
C
         DO 800 NI=1,NBF
            DO 810 NA=1,NBF
              ERR(NA,NI) = ZA(NA,NI) + ZB(NA,NI) - YA(NA,NI) - YB(NA,NI)
  810       CONTINUE
  800    CONTINUE
C
C     ---- CARRY OUT DIIS ----
C
         ETHSAV = ETHRSH
         ETHRSH = MAX(5.0D+00,ETHRSH)
         CALL CDIISU(IW,ITDIIS,ZA,ZB,ERR,WRK,A,X,
     *               IPVT,B,IODIIS,NOCP,MAXIT,MAXIT2,
     *               MAXIO,ERDIIS,NOTOPN,NBF,NBF3,1)
         ETHRSH = ETHSAV
      END IF
C
C     ---- CHECK CONVERGENCE OF SOLUTION ----
C
      DELPRV=DELTA
      DELTA=ZERO
      DO 900 NI=1,NBF
         DO 910 NA=1,NBF
            DELTA=DELTA+(ZA(NA,NI)-YA(NA,NI)+ZB(NA,NI)-YB(NA,NI))**2
  910    CONTINUE
  900 CONTINUE
      DELTA=SQRT(DELTA/NBF3)
      IF(NPRINT.NE.-5) WRITE(IW,9020) DELTA
C
C        NEXT IS CONVERGED EXIT
C
      IF(DELTA.LT.CVGTOL .AND. ITER.GE.2) THEN
         CALL RACLOS(IRAF,'DELETE')
         NOTOPN = .TRUE.
         WRITE(IW,9000) ITER
         IF(DBUG) THEN
            WRITE(IW,9030)
            CALL PRSQ(PAI1,NOA,NVIRA,NVIRA)
            WRITE(IW,9035)
            CALL PRSQ(PAI2,NOB,NVIRB,NVIRB)
         END IF
         RETURN
      END IF
C
      CALL DCOPY(NBF3,ZA,1,YA,1)
      CALL DCOPY(NBF3,ZB,1,YB,1)
C
C        CONTINUE WORKING ON CONVERGENCE
C
      IF(ITER.LT.MAXIT)  GO TO 100
C
C        THE CPHF HAS FAILED TO CONVERGE
C
      CALL RACLOS(IRAF,'DELETE')
      NOTOPN = .TRUE.
      WRITE(IW,9000) ITER
      WRITE(IW,9010) DELTA,DELPRV,CVGTOL
      CALL ABRT
      RETURN
C
 9000 FORMAT(1X,'.....',I4,' MP2 DIIS ITERATIONS CARRIED OUT .....'/)
 9010 FORMAT(1X,'CPHFU: TOO MANY ITERATIONS IN MP2 CPHF'/
     *       1X,'   FINAL ITERATION Z-MATRIX ERROR=',1P,E13.6/
     *       1X,'PREVIOUS ITERATION Z-MATRIX ERROR=',1P,E13.6/
     *       1X,' CPHF CONVERGENCE CRITERION IS',1P,E13.6/
     *       1X,'PLEASE RECONSIDER:'/
     *       5X,'IF THE SCF HAS CONVERGED TO THE RIGHT ORBITAL',
     *          ' OCCUPANCY OR NOT'/
     *       5X,'IF THIS CHEMICAL SYSTEM IS A SINGLE',
     *          ' REFERENCE WAVEFUNCTION OR NOT'/
     *       1X,'IF YOU BELIEVE THIS WAVEFUNCTION IS',
     *          ' CORRECT AND CONVERGING,'/
     *       1X,'CONSIDER INCREASING MAXITC IN $MP2.')
 9030 FORMAT(/' PAI(VIRT,OCC) ALPHA')
 9035 FORMAT(/' PAI(VIRT,OCC) BETA')
 9020 FORMAT(1X,'MAXIMUM Z-VECTOR ERROR=',1P,E13.6)
      END
C
C*MODULE MP2GRD  *DECK FNSHW2
      SUBROUTINE FNSHW2(WOCC3A,WOCC3B,PAIA,PAIB,WOCCA,WOCCB,WAIA,WAIB,
     *                  EA,EB,PMNA,PMNB,XLMNA,XLMNB,X,YA,YB,BUF,IBUF,
     *                  CA,CB,NOA,NOB,NVIRA,NVIRB,NBF,NBF2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00)
C
      LOGICAL PACK2E
C
      DIMENSION WOCC3A(NOA,NOA),PAIA(NVIRA,NOA),
     *          WOCCA(NOA,NOA),WAIA(NVIRA,NOA),CA(NBF,NBF)
      DIMENSION WOCC3B(NOB,NOB),PAIB(NVIRB,NOB),
     *          WOCCB(NOB,NOB),WAIB(NVIRB,NOB),CB(NBF,NBF)
      DIMENSION XLMNA(NBF,NBF),PMNA(NBF2),EA(*),YA(NBF,NBF),
     *          XLMNB(NBF,NBF),PMNB(NBF2),EB(*),YB(NBF,NBF),
     *          X(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
C
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
C
C     ---- ADD FINAL TERMS TO THE ENERGY WEIGHTED DENSITY ----
C
      CALL VCLR(WOCC3A,1,NOA*NOA)
      CALL VCLR(WOCC3B,1,NOB*NOB)
C
C     ---- FORM FOCK-LIKE MATRICES YA(MU,NU) AND YB(MU,NU) ----
C
      CALL FLMATU(PMNA,YA,XLMNA,PMNB,YB,XLMNB,BUF,IBUF,NBF)
C
C     ---- TRANSFORM THE MATRICES YA AND YB TO MO BASIS ----
C
      CALL MRARBR(YA,NBF,NBF,NBF,CA,NBF,NOA,X,NBF)
      CALL DGEMM('T','N',NOA,NOA,NBF,ONE,CA,NBF,X,NBF,
     *           ONE,WOCC3A,NOA)
      CALL MRARBR(YB,NBF,NBF,NBF,CB,NBF,NOB,X,NBF)
      CALL DGEMM('T','N',NOB,NOB,NBF,ONE,CB,NBF,X,NBF,
     *           ONE,WOCC3B,NOB)
C
C     ---- COMPLETE THE ENERGY WEIGHTED DENSITY IN THE MO BASIS ----
C
      DO 50 J=1,NOA
         DO 60 I=1,NOA
            WOCCA(I,J)=WOCCA(I,J) - WOCC3A(I,J)
   60    CONTINUE
   50 CONTINUE
      DO 70 J=1,NOB
         DO 80 I=1,NOB
            WOCCB(I,J)=WOCCB(I,J) - WOCC3B(I,J)
   80    CONTINUE
   70 CONTINUE
C
C     ---- ADD ENERGY TERM TO WAI ----
C
      DO 125 I=1,NOA
         DO 120 J=1,NVIRA
            WAIA(J,I) = WAIA(J,I) - EA(I)*PAIA(J,I)
  120    CONTINUE
  125 CONTINUE
      DO 135 I=1,NOB
         DO 130 J=1,NVIRB
            WAIB(J,I) = WAIB(J,I) - EB(I)*PAIB(J,I)
  130    CONTINUE
  135 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UMPGTR
      SUBROUTINE UMPGTR(PIJ,PAB,PAI,C,NBF,NOC,NVIR,
     *           X,Y,PMN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (HALF=0.5D+00,ONE=1.0D+00,TWO=2.0D+00)
C
      DIMENSION PIJ(NOC,NOC),PAB(NVIR,NVIR),PAI(NVIR,NOC),
     *          C(NBF,*),X(NBF,*),Y(NBF,*),PMN(*)
C
      CALL MRARTR(PIJ,NOC,NOC,NOC,C,NBF,NBF,X,NBF)
      CALL MRARBR(C,NBF,NBF,NOC,X,NBF,NBF,Y,NBF)
C
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,ONE,
     *     Y,NBF)
C
      CALL MRARTR(PAI,NVIR,NVIR,NOC,C,NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,TWO,C(1,NOC+1),NBF,X,NBF,ONE,
     *     Y,NBF)
C
      MUNU=0
      DO 15 MU=1,NBF
         DO 10 NU=1,MU
            MUNU=MUNU+1
            PMN(MUNU) =HALF*(Y(MU,NU) + Y(NU,MU))
   10    CONTINUE
   15 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK UMPW2
      SUBROUTINE UMPW2(V,E,D1,D2,NOC,N1,KRECX7,KRECX8,KRECX9,LRHO)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C     COMMON /DAFILE/ IDAF,NAV,IODA(255)
C
      DIMENSION D1(*),D2(*),V(N1,*),E(*)
C
      N2=(N1*N1+N1)/2
C
C     -----  D2 IS THE MP2 ENERGY-WEIGHTED DENSITY  -----
C
      OCC=ONE
      IJ=0
      DO 25 I=1, N1
         DO 20 J=1, I
            IJ=IJ+1
            DUM=ZERO
            DO 30 K=1, NOC
               DUM=DUM-E(K)*V(I,K)*V(J,K)
   30       CONTINUE
            D2(IJ)=D2(IJ)+DUM*OCC
   20    CONTINUE
   25 CONTINUE
C
C     -----  D2 IS NOW HF+MP2 ENERGY-WEIGHTED DENSITY -----
C
      CALL DAWRIT(IDAF,IODA,D2,N2,KRECX9,0)
C
C     -----  NOW FORM TOTAL HF+MP2 DENSITY  -----
C
      CALL DAREAD(IDAF,IODA,D1,N2,LRHO,0)
      CALL DAREAD(IDAF,IODA,D2,N2,KRECX7,0)
      DO 10 N=1, N2
         D2(N)=D1(N) + D2(N)
   10 CONTINUE
C
      CALL DAWRIT(IDAF,IODA,D2,N2,LRHO,0)
C
C     ----- SAVE PHF ON KRECX8 -----
C
      CALL DAWRIT(IDAF,IODA,D1,N2,KRECX8,0)
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK NONSEP
      SUBROUTINE NONSEP(AIBJ,AIBJ2,EA,EB,CA,CB,X,Y,XX,XIN,NBF,NBF2,NOA,
     *                  NOB,NVIRA,NVIRB,NFTAA,NFTAB,NFT16,TOL,SOME,
     *                  NACORE,NBCORE,NOCA,NOCB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL SOME
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /IJPAIR/ IA(MXAO)
C
      DIMENSION AIBJ(NVIRA,NVIRA,NOCA),AIBJ2(NVIRB,NVIRA,NOCB),
     *          CA(NBF,*),CB(NBF,*),X(NBF,*),Y(NBF,*),
     *          XIN(NBF2,NBF),XX(*),EA(*),EB(*)
C
C         GENERATE NON-SEPARABLE CONTRIBUTIONS TO MP2 2ND ORDER DENSITY
C         PROGRAM ADAPTED FROM MAMNLS
C
      CALL SEQREW(NFTAA)
      CALL SEQREW(NFTAB)
      NO = 0
C
      DO 1000 MI=1,NOCA
C
C     ----- READ ALL (AI|BJ) INTEGRALS FOR GIVEN -I- -----
C
      DO 115 MA=1,NVIRA
         CALL TRFRD(NFTAA,XX,NX,X,NIX,NBF2)
         DO 110 MB=1,NVIRA
            MBJ=IA(MB+NOA)+NACORE
            DO 105 MJ=1,NOCA
               AIBJ(MB,MA,MJ) = XX(MBJ+MJ)
  105       CONTINUE
  110    CONTINUE
  115 CONTINUE
C
      CALL VCLR(XIN,1,NBF*NBF2)
C
C     ----- -X(A,B)- IS -GAMMA(A,I,B,J)- FOR GIVEN -I,J- -----
C
      DO 150 MJ=1,NOCA
         DO 135 MB=1,NVIRA
            DO 130 MA=1,NVIRA
               DIJAB = EA(MI+NACORE) + EA(MJ+NACORE)
     *                           - EA(MA+NOA) - EA(MB+NOA)
               X(MA,MB) = (AIBJ(MB,MA,MJ)-
     *                    AIBJ(MA,MB,MJ))/DIJAB
  130       CONTINUE
  135    CONTINUE
C         CALL PRSQ(X,NVIRA,NVIRA,NVIRA)
C         CALL PRSQ(CA,NBF,NBF,NBF)
C
C     ----- TRANSFORM -A,B- TO -MU,LAMBDA- -----
C
         CALL MRARTR(X,NBF,NVIRA,NVIRA,CA(1,NOA+1),NBF,NBF,Y,NBF)
         CALL MRARBR(CA(1,NOA+1),NBF,NBF,NVIRA,Y,NBF,NBF,X,NBF)
C
C     ----- TRANSFORM -J- TO -SIGMA- AND SYMMETRIZE -SIGMA,LAMBDA -----
C     ----- -XIN(MU,LASI)- IS -GAMMA(MU,I,LASI)- FOR GIVEN -I-    -----
C
         DO 148 NMU=1, NBF
            LASI = 0
            DO 143 NLA=1, NBF
               DO 140 NSI=1, NLA
                  LASI = LASI+1
                  XIN(LASI,NMU)=XIN(LASI,NMU)
     *                          + CA(NSI,MJ+NACORE)*X(NMU,NLA)
     *                          + CA(NLA,MJ+NACORE)*X(NMU,NSI)
  140          CONTINUE
  143       CONTINUE
  148   CONTINUE
  150 CONTINUE
C
C PREPARE CROSS TERM IN SIMILAR FASHION
C
      DO 215 MA=1,NVIRA
         CALL TRFRD(NFTAB,XX,NX,X,NIX,NBF2)
         DO 210 MB=1,NVIRB
            MBJ=IA(MB+NOB)+NBCORE
            DO 205 MJ=1,NOCB
               AIBJ2(MB,MA,MJ) = XX(MBJ+MJ)
  205       CONTINUE
  210    CONTINUE
  215 CONTINUE
C
      DO 250 MJ=1,NOCB
         DO 230 MA=1,NVIRA
            DO 225 MB=1,NVIRB
               DIJAB = EA(MI+NACORE) + EB(MJ+NBCORE)
     *                           - EA(MA+NOA) - EB(MB+NOB)
               X(MA,MB) = AIBJ2(MB,MA,MJ)/DIJAB
  225       CONTINUE
  230    CONTINUE
C         CALL PRSQ(X,NVIRB,NVIRA,NVIRA)
C         CALL PRSQ(CA,NBF,NBF,NBF)
C         CALL PRSQ(CB,NBF,NBF,NBF)
C
C     ----- TRANSFORM -A,B- TO -MU,LAMBDA- -----
C
         CALL MRARTR(X,NBF,NVIRA,NVIRB,CB(1,NOB+1),NBF,NBF,Y,NBF)
         CALL MRARBR(CA(1,NOA+1),NBF,NBF,NVIRA,Y,NBF,NBF,X,NBF)
C
C     ----- TRANSFORM -J- TO -SIGMA- AND SYMMETRIZE -SIGMA,LAMBDA -----
C     ----- -XIN(MU,LASI)- IS -GAMMA(MU,I,LASI)- FOR GIVEN -I-    -----
C
         DO 248 NMU=1, NBF
            LASI = 0
            DO 243 NLA=1, NBF
               DO 240 NSI=1, NLA
                  LASI = LASI+1
                  XIN(LASI,NMU)=XIN(LASI,NMU)
     *                          + CB(NSI,MJ+NBCORE)*X(NMU,NLA)
     *                          + CB(NLA,MJ+NBCORE)*X(NMU,NSI)
  240          CONTINUE
  243       CONTINUE
  248   CONTINUE
  250 CONTINUE
C
      DO 160 NMU=1, NBF
         CALL TRFWTM(NFT16,XIN(1,NMU),NX,X,NIX,NBF2,1,TOL)
         NO = NO+NX
  160 CONTINUE
C
 1000 CONTINUE
C
      IF(SOME) WRITE(IW,9999) NO,NFT16
      RETURN
C
 9999 FORMAT(1X,I12,' -2DM- ELEMENTS WRITTEN ON FILE',I3)
      END
C
C*MODULE MP2GRD  *DECK JKDUMP
      SUBROUTINE JKDUMP(PHFA,PMP2A,CA,PHFB,PMP2B,CB,XINA,XINB,
     *                  NBF3,NBF2,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2DM2/ NOC1,NOC2,ISTEP,NOC3,NOC4
C
      DIMENSION PHFA(NBF2),PMP2A(NBF2),XINA(NBF2,NBF,*),CA(NBF3)
      DIMENSION PHFB(NBF2),PMP2B(NBF2),XINB(NBF2,NBF,*),CB(NBF3)
C
      DATA KRECA7,KRECA8,KRECB7,KRECB8 /417,418,427,428/
      DATA NFT16,NFT15/16,15/
C
C        READ PARTIALLY TRANSFORMED -DM2- FROM DISK FOR GIVEN I RANGE
C        NOTE: PHFA ONLY THERE AS A DUMMY
C
      IF(ISTEP.EQ.1) THEN
         CALL SEQREW(NFT16)
         CALL SEQREW(NFT15)
      END IF
      DO 15 MI=1,NOC2-NOC1+1
         DO 10 NMU=1,NBF
            CALL TRFRD(NFT16,XINA(1,NMU,MI),NX,PHFA,NIX,NBF2)
   10    CONTINUE
   15 CONTINUE
      DO 25 MI=1,NOC4-NOC3+1
         DO 20 NMU=1,NBF
            CALL TRFRD(NFT15,XINB(1,NMU,MI),NX,PHFA,NIX,NBF2)
   20    CONTINUE
   25 CONTINUE
C
C        READ SCF AND MP2 FIRST ORDER DENSITY MATRICES
C
      IF(ISTEP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,PHFA ,NBF2,KRECA8,0)
         CALL DAREAD(IDAF,IODA,PMP2A,NBF2,KRECA7,0)
         CALL DAREAD(IDAF,IODA,CA   ,NBF3,   15,0)
         CALL DAREAD(IDAF,IODA,PHFB ,NBF2,KRECB8,0)
         CALL DAREAD(IDAF,IODA,PMP2B,NBF2,KRECB7,0)
         CALL DAREAD(IDAF,IODA,CB   ,NBF3,    19,0)
      END IF
      RETURN
      END
