C 11 APR 08 - AA  - COSMETIC REFORMATTING OF COMMON BLOCKS
C 22 DEC 06 - ST,NK,MC - ADD LC EXCHANGE ARGUMENTS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  5 FEB 05 - MWS - CHANGE COMMON NAME GAMMAF TO BE FMTTBL INSTEAD
C 10 NOV 04 - KI  - MCDAV,INTJ4-6,INTK5,INTK6: FASTER SP ROTATED AXIS
C  7 SEP 04 - GDF - INTJ1-INTJ6: USE FULLY ACCURATE FM(T) INTERPOLATION
C  4 NOV 03 - KI  - SP0S1S: DO AXES-1 INTEGRALS BY MCMURCHIE-DAVIDSON,
C                   INCORPORATE ROTATION TO AXES2; GENR70,TQ0S1S: TWEAKS
C                   INTJ1 TO INTJ6: NEW AUXILIARIES FOR M-D INTEGRATION
C  7 AUG 02 - JMS - ALTERED COMMENTS TO REFLECT GRADIENT CHANGES
C 16 NOV 01 - JMS - GAMGEN HERE, POPLE PACKAGE INITIALIZATIONS CHANGED
C  8 OCT 01 - JMS - ADDITIONAL REVISIONS
C  1 AUG 01 - JMS - COMPLETE REWRITE OF ROTATED AXIS INTEGRAL CODE.
C                   HISTORY BELOW IS KEPT MAINLY FOR HISTORICAL REASONS
C 13 JUN 01 - JMS - PI/4 VALUE CHANGED FROM INACCURATE 0.7853981625
C 13 JUN 96 - MWS - REMOVE FTNCHEK WARNINGS
C 10 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 15 JUL 94 - BMB - REDIMENSIONED TABLE TO 7200 TOTAL ELEMENTS
C 25 SEP 91 - MWS - COMMON GOUT CONTAINS ONLY POPLE METHOD INTEGRALS
C  9 AUG 91 - TLW - REDIMENSIONED GOUT FROM 3888 TO 151875
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 20 MAY 88 - MWS - VECTORIZE TQ1111, PASS GOUT TO TQ0111 BY CALL
C  9 MAY 88 - MWS - USE PARAMETER DIMENSIONS FOR COMMON BLOCKS
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTATNS
C  4 APR 86 - MWS - LOGICAL OUT IN SP0000
C 22 OCT 85 - STE - USE GENERIC EXP,SQRT; ROT2: IMPROVE LOGIC;
C                   SP0000: USE /CHLT/TOL IN EXPONENT SIZE CHECK
C                   REMOVE DATA DEPENENCIES IN R30011,R30101,R30111
C                   R31111,TQ1111
C 10 JAN 84 - STE - ROLL ROT2 CODE BACK INTO LOOPS
C 16 DEC 83 - STE - ROLL TQ1111 CODE BACK INTO LOOPS
C  2 OCT 82 - MWS,NDSU - CONVERT TO IBM
C
C*MODULE INT2B   *DECK GENR70
      SUBROUTINE GENR70(IEXCH,SCHSKP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SCHSKP
C
C     GAUSSIAN TWO ELECTRON INTEGRAL PACKAGE
C
C     SCHSKP (SCHWARZ SKIP) MEANS ONLY IB INDEXING, AND ZEROING IS DONE.
C
C     GENR70 EVALUATES THOSE INTEGRALS INVOLVING ONLY S AND P FUNCTIONS.
C     MAIN LOOP OVER SHELLS ... SEE EXPLANATION IN ROUTINE GINPUT ...
C     LINK 301
C     ACCEPTS NUMBERS OF FOUR SHELLS: ISHELL  JSHELL  KSHELL  LSHELL
C     FINDS THEIR ANGULAR QUANTUM NUMBERS AND, BASED ON THIS
C     ORDERS SHELLS IN A STANDARD MANNER: INEW  JNEW  KNEW  LNEW
C
C     ONLY POSSIBILITIES ALLOWED FOR ANGULAR QUANTUM NUMBERS ARE THEN
C     0000  0001  0011  0101  0111  1111
C     DETERMINES TYPE OF INTEGRAL SET BASED ON THE ABOVE NUMBERS
C     CALLS THE FOLLOWING ROUTINES IN THE ORDER GIVEN (THESE ROUTINES
C     WERE LIFTED FROM GAUSSIAN76):
C
C     SINFO, TO PRESET INTEGRAL ACCURACY LIMITS
C
C     SGEOM, OBTAINS GEOMETRICAL INFORMATION ABOUT THE FOUR CENTERS
C     AND FINDS TWO SETS OF LOCAL AXES FOR CENTERS:
C     (A AND B) P SET AND (C AND D) Q SET
C
C     PINF, OBTAINS INFORMATION ABOUT GAUSSIAN FUNCTIONS CONNECTED
C     WITH THE P SET OF AXES
C
C     AT THIS POINT GENR70 OBTAINS INFORMATION ABOUT THE GAUSSIAN
C     FUNCTIONS CONNECTED WITH THE Q SET OF AXES
C
C     SP0000 TO SP1111, OBTAIN UP TO 88 INTEGRALS REFERRED
C     TO AXES A B AND Q
C
C     ROT2, ROTATES THESE INTEGRALS TO UP TO 160 INTEGRALS ON A B AND Q
C
C     TQ0011 TO TQ1111, TRANSLATE THESE INTEGRALS ON A B AND Q
C     TO UP TO 256 INTEGRALS ON A B C AND D
C
C     R30001 TO R31111, ROTATE UP TO 256 INTEGRALS ON A B C AND D
C     TO THE SAME NUMBER REFERRED TO THE FIXED SPACE AXES
C
C          THE POPLE/HEHRE ROTATED AXIS CODE WAS EXTENSIVELY REWRITTEN
C          IN JULY 2001 BY JOSE M. SIERRA OF SYNSTAR, SA, MADRID SPAIN.
C     THE CHANGES INCLUDE USING DO LOOPS INSTEAD OF UNROLLED CODE, AND
C     ROUTINES SINFO, SGEOM AND PINF HAVE BEEN IN-LINED INTO GENR70
C     ROUTINE  R30S1S ENCOMPASSES SUBS FROM R30001 TO R31111
C     ROUTINE  SP0S1S ENCOMPASSES SUBS FROM SP0000 TO SP1111
C     ROUTINE  TQ0S1S ENCOMPASSES SUBS FROM TQ0011 TO TQ1111
C
      PARAMETER (MXGTOT=20000, MXSH=5000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      COMMON /B     / CO(MXSH,3)
      COMMON /FLIPS / IB(4,3)
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /GOUT  / GPOPLE(256*3),NORG
      COMMON /INTAC2/ EI1,EI2,CUX
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     2                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     3                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     4                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /POPOUT/ LPOPI,LPOPJ,LPOPK,LPOPL
      COMMON /SHLG70/ ISHELL,JSHELL,KSHELL,LSHELL,INEW,JNEW,KNEW,LNEW
      COMMON /SHLLFO/ NGA,LA,EXA(MXGSH),CSA(MXGSH),CPA(MXGSH),
     2                NGB,LB,EXB(MXGSH),CSB(MXGSH),CPB(MXGSH),
     3                NGC,LC,EXC(MXGSH),CSC(MXGSH),CPC(MXGSH),
     4                NGD,LD,EXD(MXGSH),CSD(MXGSH),CPD(MXGSH)
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     2                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     3                NIJ,IJ,KL,IJKL
C
C     ==================================================================
C
CJMS  LABELLED COMMON JMSGYH DEFINED FOR COMPUTATIONAL EFFICIENCY.
CJMS  IT IS ONLY USED IN THIS MODULE INT2B.
C
      COMMON /JMSGYH/ SQ(0:1,0:1)
      COMMON /KI3   / R00(25),R01(120),R02(156),R03(80),R04(15)
      COMMON /KI2   / ACY,ACY2,AQX,AQX2,AQXY,Y03,Y04
C
      DIMENSION P12(3,3),P34(3,3),P(3,3),T(3)
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00, PT5=0.5D+00)
      PARAMETER (PT7=0.7D+00, PT9=0.9D+00)
C        ACYCUT CHANGED FROM A VALUE OF 1E-4 IN 2004,
C        FOR ACCURACY WHEN USING DIFFUSE EXPONENTS.
      PARAMETER (ACYCUT=1.0D-10, TENM12=1.0D-12)
C
CJMS  PARAMETER (PIDIV4=0.785398163397448279D+00)
CC    PARAMETER (PIDIV4=0.785398163397448309615660845819D+00)
CC    PARAMETER (PIDIV4=0.785398163397448D+00)
      PARAMETER (PITO52=0.349868366552497D+02)
C
CC    PIDIV4= ATAN(ONE)
CC    PI4= PIDIV4
CC    PI = PI4+PI4+PI4+PI4
CC    PITO52=(PI +PI )*PI * SQRT(PI )
      PIF= PITO52
C
      LAT= KTYPE(ISHELL)-1
      LBT= KTYPE(JSHELL)-1
      LCT= KTYPE(KSHELL)-1
      LDT= KTYPE(LSHELL)-1
      ITYPE= 1+LDT+2*(LCT+2*(LBT+2*LAT))
C
      IF(ITYPE.EQ. 1 .OR. ITYPE.EQ. 2 .OR. ITYPE.EQ. 4 .OR.
     2   ITYPE.EQ. 6 .OR. ITYPE.EQ. 8 .OR. ITYPE.EQ.16) THEN
C
C TYPES 0000,0001,0011,0101,0111,1111 ARE UNALTERED
C
         INEW= ISHELL
         JNEW= JSHELL
         KNEW= KSHELL
         LNEW= LSHELL
         LA= LAT
         LB= LBT
         LC= LCT
         LD= LDT
         IB(1,IEXCH)= 1
         IB(2,IEXCH)= 2
         IB(3,IEXCH)= 3
         IB(4,IEXCH)= 4
C
         LPOPI = 64
         LPOPJ = 16
         LPOPK =  4
         LPOPL =  1
      ELSEIF(ITYPE.EQ.10 .OR. ITYPE.EQ.12) THEN
C
C TYPES 1001,1011 HAVE IJ SWITCHED
C
         INEW= JSHELL
         JNEW= ISHELL
         KNEW= KSHELL
         LNEW= LSHELL
         LA= LBT
         LB= LAT
         LC= LCT
         LD= LDT
         IB(1,IEXCH)= 2
         IB(2,IEXCH)= 1
         IB(3,IEXCH)= 3
         IB(4,IEXCH)= 4
C
         LPOPI = 16
         LPOPJ = 64
         LPOPK =  4
         LPOPL =  1
      ELSEIF(ITYPE.EQ. 3 .OR. ITYPE.EQ. 7) THEN
C
C TYPES 0010,0110 HAVE KL SWITCHED
C
         INEW= ISHELL
         JNEW= JSHELL
         KNEW= LSHELL
         LNEW= KSHELL
         LA= LAT
         LB= LBT
         LC= LDT
         LD= LCT
         IB(1,IEXCH)= 1
         IB(2,IEXCH)= 2
         IB(3,IEXCH)= 4
         IB(4,IEXCH)= 3
C
         LPOPI = 64
         LPOPJ = 16
         LPOPK =  1
         LPOPL =  4
      ELSEIF(ITYPE.EQ. 5 .OR. ITYPE.EQ.13 .OR. ITYPE.EQ.14) THEN
C
C TYPES 0100,1100,1101 HAVE PAIRS IJ AND KL SWITCHED
C
         INEW= KSHELL
         JNEW= LSHELL
         KNEW= ISHELL
         LNEW= JSHELL
         LA= LCT
         LB= LDT
         LC= LAT
         LD= LBT
         IB(1,IEXCH)= 3
         IB(2,IEXCH)= 4
         IB(3,IEXCH)= 1
         IB(4,IEXCH)= 2
C
         LPOPI =  4
         LPOPJ =  1
         LPOPK = 64
         LPOPL = 16
      ELSEIF(ITYPE.EQ.11) THEN
C
C TYPE 1010 HAS IJ SWITCHED AND KL SWITCHED
C
         INEW= JSHELL
         JNEW= ISHELL
         KNEW= LSHELL
         LNEW= KSHELL
         LA= LBT
         LB= LAT
         LC= LDT
         LD= LCT
         IB(1,IEXCH)= 2
         IB(2,IEXCH)= 1
         IB(3,IEXCH)= 4
         IB(4,IEXCH)= 3
C
         LPOPI = 16
         LPOPJ = 64
         LPOPK =  1
         LPOPL =  4
      ELSEIF(ITYPE.EQ. 9) THEN
C
C TYPE 1000  HAS PAIRS IJ AND KL SWITCHED FOLLOWED BY KL SWITCH
C
         INEW= KSHELL
         JNEW= LSHELL
         KNEW= JSHELL
         LNEW= ISHELL
         LA= LCT
         LB= LDT
         LC= LBT
         LD= LAT
         IB(1,IEXCH)= 4
         IB(2,IEXCH)= 3
         IB(3,IEXCH)= 1
         IB(4,IEXCH)= 2
C
         LPOPI =  1
         LPOPJ =  4
         LPOPK = 64
         LPOPL = 16
      ELSEIF(ITYPE.EQ.15) THEN
C
C TYPE 1110 HAS PAIRS IJ AND KL SWITCHED FOLLOWED BY IJ SWITCH
C
         INEW= LSHELL
         JNEW= KSHELL
         KNEW= ISHELL
         LNEW= JSHELL
         LA= LDT
         LB= LCT
         LC= LAT
         LD= LBT
         IB(1,IEXCH)= 3
         IB(2,IEXCH)= 4
         IB(3,IEXCH)= 2
         IB(4,IEXCH)= 1
C
         LPOPI =  4
         LPOPJ =  1
         LPOPK = 16
         LPOPL = 64
      ENDIF
C
C ONLY 6 STANDARD TYPES REMAIN. 0000,0001,0011,0101,0111,1111
C SPECIFY THESE BY JTYPE
C
      IF(ITYPE.EQ. 1) THEN
         JTYPE= 1
         R00(1)= ZER
      ELSEIF(ITYPE.EQ. 2 .OR. ITYPE.EQ. 3 .OR.
     2       ITYPE.EQ. 5 .OR. ITYPE.EQ. 9) THEN
         JTYPE= 2
         DO 310 I=1,2
            R00(I)=ZER
  310    CONTINUE
         DO 320 I=1,3
            R01(I)=ZER
  320    CONTINUE
      ELSEIF(ITYPE.EQ. 4 .OR. ITYPE.EQ.13) THEN
         JTYPE= 3
         DO 330 I=1,5
            R00(I)=ZER
  330    CONTINUE
         DO 340 I=1,12
            R01(I)=ZER
  340    CONTINUE
         DO 350 I=1,6
            R02(I)=ZER
  350    CONTINUE
      ELSEIF(ITYPE.EQ. 6 .OR. ITYPE.EQ. 7 .OR.
     2       ITYPE.EQ.10 .OR. ITYPE.EQ.11) THEN
         JTYPE= 4
         DO 360 I=1,4
            R00(I)=ZER
  360    CONTINUE
         DO 370 I=1,12
            R01(I)=ZER
  370    CONTINUE
         DO 380 I=1,6
            R02(I)=ZER
  380    CONTINUE
      ELSEIF(ITYPE.EQ. 8 .OR. ITYPE.EQ.12 .OR.
     2       ITYPE.EQ.14 .OR. ITYPE.EQ.15) THEN
         JTYPE= 5
         DO 390 I=1,10
            R00(I)=ZER
  390    CONTINUE
         DO 400 I=1,39
            R01(I)=ZER
  400    CONTINUE
         DO 410 I=1,36
            R02(I)=ZER
  410    CONTINUE
         DO 420 I=1,10
            R03(I)=ZER
  420    CONTINUE
      ELSEIF(ITYPE.EQ.16) THEN
         JTYPE= 6
         DO 430 I=1,25
            R00(I)=ZER
  430    CONTINUE
         DO 440 I=1,120
            R01(I)=ZER
  440    CONTINUE
         DO 450 I=1,156
            R02(I)=ZER
  450    CONTINUE
         DO 460 I=1,80
            R03(I)=ZER
  460    CONTINUE
         DO 470 I=1,15
            R04(I)=ZER
  470    CONTINUE
      ENDIF
C
C EMPTY COMMON GOUT
C
C---  NGOUT=256
C---  CALL VCLR(GPOPLE(1+NORG),1,NGOUT)
C
      IF(SCHSKP) GO TO 999
C
C OBTAIN INFORMATION ABOUT SHELLS: INEW, KNEW, JNEW, LNEW
C NUMBER OF GAUSSIANS GO INTO NGA,... IN COMMON SHLLFO
C SHELL ANGULAR QUANTUM NUMBERS LA,... GO INTO COMMON SHLLFO
C GAUSSIAN EXPONENTS GO INTO ARRAYS EXA,EXB,EXC,EXD IN COMMON SHLLFO
C GAUSSIAN COEFFICIENTS GO INTO ARRAYS CSA,CPA,... IN COMMON SHLLFO
C
C NUMBERS OF GAUSSIAN FUNCTIONS IN SHELLS INEW JNEW KNEW AND LNEW
C
      NGA= KNG(INEW)
      NGB= KNG(JNEW)
      NGC= KNG(KNEW)
      NGD= KNG(LNEW)
C
C STARTING LOCATIONS OF SHELLS INEW JNEW KNEW AND LNEW IN LIST
C OF GAUSSIAN FUNCTIONS
C
      I= KSTART(INEW)-1
      J= KSTART(JNEW)-1
      K= KSTART(KNEW)-1
      L= KSTART(LNEW)-1
C
C LOOP OVER GAUSSIANS IN EACH SHELL
C FIRST SHELL INEW
C
      DO 110 NI=1,NGA
         N=I+NI
C
C THE MAXIMUM COEFFICIENT ASSOCIATED WITH SHELL IS
C USED TO DETERMINE IF ANY OF THE INTEGRALS ASSOCIATED WITH A SET
C OF SHELLS IS LARGE ENOUGH TO WARRANT EVALUATION OF THE ENTIRE SET
C
         CMAXA(NI)= CMAX(N)
C
C GAUSSIAN EXPONENTS
C
         EXA(NI)= EX(N)
C
C S COEFFICIENTS
C
         CSA(NI)= CS(N)
C
C P COEFFICIENTS
C
         CPA(NI)= CP(N)
  110 CONTINUE
C
C REPEAT PROCEDURE FOR SHELLS JNEW KNEW AND LNEW
C
      DO 120 NJ=1,NGB
         N=J+NJ
         CMAXB(NJ)= CMAX(N)
         EXB(NJ)= EX(N)
         CSB(NJ)= CS(N)
         CPB(NJ)= CP(N)
  120 CONTINUE
      DO 130 NK=1,NGC
         N=K+NK
         CMAXC(NK)=CMAX(N)*QQ4
         EXC(NK)=EX(N)
         CSC(NK)=CS(N)*QQ4
         CPC(NK)=CP(N)*QQ4
  130 CONTINUE
      DO 140 NL=1,NGD
         N=L+NL
         CMAXD(NL)=CMAX(N)
         EXD(NL)=EX(N)
         CSD(NL)=CS(N)
         CPD(NL)=CP(N)
  140 CONTINUE
C
C     MAB= LA+LB-1
C     MCD= LC+LD-1
      NGANGB=NGA*NGB
C
C COORDINATES OF ATOMS ASSOCIATED WITH SHELLS INEW JNEW KNEW AND LNEW
C
      R12= ZER
      R34= ZER
      DO 150 N=1,3
         P12(N,1)= CO(INEW,N)
         P12(N,2)= CO(JNEW,N)
         P12(N,3)= P12(N,2)-P12(N,1)
      R12= R12+P12(N,3)*P12(N,3)
         P34(N,1)= CO(KNEW,N)
         P34(N,2)= CO(LNEW,N)
         P34(N,3)= P34(N,2)-P34(N,1)
  150 R34= R34+P34(N,3)*P34(N,3)
C
C FIND DIRECTION COSINES OF PENULTIMATE AXES FROM COORDINATES OF AB
C P(1,1),P(1,2),... ARE DIRECTION COSINES OF AXES AT P.  Z-AXIS ALONG AB
C T(1),T(2),T(3)... ARE DIRECTION COSINES OF AXES AT Q.  Z-AXIS ALONG CD
C
C FIND DIRECTION COSINES OF AB AND CD. THESE ARE LOCAL Z-AXES.
C IF INDETERMINATE TAKE ALONG SPACE Z-AXIS
C
      P(1,3)= ZER
      P(2,3)= ZER
      P(3,3)= ONE
      RAB= ZER
      IF(R12.NE.ZER) THEN
         RAB= SQRT(R12)
         TMP= ONE/RAB
         P(1,3)= P12(1,3)*TMP
         P(2,3)= P12(2,3)*TMP
         P(3,3)= P12(3,3)*TMP
      ENDIF
C
      T(1)= ZER
      T(2)= ZER
      T(3)= ONE
      RCD= ZER
      IF(R34.NE.ZER) THEN
         RCD= SQRT(R34)
         TMP= ONE/RCD
         T(1)= P34(1,3)*TMP
         T(2)= P34(2,3)*TMP
         T(3)= P34(3,3)*TMP
      ENDIF
C
C FIND LOCAL Y-AXIS AS COMMON PERPENDICULAR TO AB AND CD
C IF INDETERMINATE TAKE PERPENDICULAR TO AB AND SPACE Z-AXIS
C IF STILL INDETERMINATE TAKE PERPENDICULAR TO AB AND SPACE X-AXIS
C
      COSG= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
      COSG= MIN( ONE,COSG)
      COSG= MAX(-ONE,COSG)
C     SING= SQRT(ONE-COSG*COSG)
C
C MODIFIED ROTATION TESTING.
C THIS FIX CURES THE SMALL ANGLE PROBLEM.
C
      P(1,2)= T(3)*P(2,3)-T(2)*P(3,3)
      P(2,2)= T(1)*P(3,3)-T(3)*P(1,3)
      P(3,2)= T(2)*P(1,3)-T(1)*P(2,3)
      IF( ABS(COSG).GT.PT9) THEN
         SING= SQRT(P(1,2)*P(1,2)+P(2,2)*P(2,2)+P(3,2)*P(3,2))
      ELSE
         SING= SQRT(ONE-COSG*COSG)
      ENDIF
      IF( ABS(COSG).LE.PT9 .OR. SING.GE.TENM12) THEN
         TMP= ONE/SING
         P(1,2)= P(1,2)*TMP
         P(2,2)= P(2,2)*TMP
         P(3,2)= P(3,2)*TMP
      ELSE
         I=3
         IF( ABS(P(1,3)).LE.PT7) I=1
         TMP = P(I,3)*P(I,3)
         TMP = MIN( ONE,TMP)
         TMP = SQRT(ONE-TMP)
         IF(TMP.NE.ZER) TMP= ONE/TMP
         IF( ABS(P(1,3)).LE.PT7) THEN
            P(1,2)= ZER
            P(2,2)= P(3,3)*TMP
            P(3,2)=-P(2,3)*TMP
         ELSE
            P(1,2)= P(2,3)*TMP
            P(2,2)=-P(1,3)*TMP
            P(3,2)= ZER
         ENDIF
      ENDIF
C
C FIND DIRECTION COSINES OF LOCAL X-AXES
C
      P(1,1)= P(2,2)*P(3,3)-P(3,2)*P(2,3)
      P(2,1)= P(3,2)*P(1,3)-P(1,2)*P(3,3)
      P(3,1)= P(1,2)*P(2,3)-P(2,2)*P(1,3)
C
C FIND COORDINATES OF C RELATIVE TO LOCAL AXES AT A
C
      T(1)= P34(1,1)-P12(1,1)
      T(2)= P34(2,1)-P12(2,1)
      T(3)= P34(3,1)-P12(3,1)
      ACX = T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
      ACY = T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
      ACZ = T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
C
C SET ACY= 0  IF CLOSE
C
      IF( ABS(ACY).LE.ACYCUT) THEN
         ACY =ZER
         ACY2=ZER
      ELSE
         ACY2=ACY*ACY
      ENDIF
C
C DIRECTION COSINES OF CD LOCAL AXES WITH RESPECT TO AB LOCAL AXES
C (COSG,0,-SING)  (0,1,0)  (SING,0,COSG)
C
C PRELIMINARY P LOOP
C
C FILL GEOMPQ WITH INFORMATION ABOUT P IN PRELIMINARY P-LOOP
C
      JI= 1
      DO 170 I=1,NGA
         X01= EXA(I)
         DO 170 J=1,NGB
            X02= EXB(J)
            X12= X01+X02
            X21= ONE/X12
            Y01= X01*X21
            Y02= ONE-Y01
            Y12= Y01*X02
            TX12(JI)= X12
            TX21(JI)= X21*PT5
            TY02(JI)= Y02*RAB
CJMS        TY01(JI)=-Y01*RAB
            TY01(JI)= TY02(JI)-RAB
C           IF(JTYPE.EQ.4 .OR. JTYPE.EQ.5) TY01(JI)= TY01(JI)*X12
            R12Y12= R12*Y12
            IF(R12Y12.GT.CUX) THEN
               ISMLP(JI)=2
               GO TO 160
            ENDIF
            E12= X21* EXP(-R12Y12)
            TST= E12*CMAXA(I)*CMAXB(J)
            ISMLP(JI)=0
            IF(TST.LE.EI1) ISMLP(JI)=1
            IF(TST.LE.EI2) ISMLP(JI)=2
            E12= PIF*E12
C
C FOR TYPES 0000,0001,0011 ONLY D00P NEEDED
C
            D00P(JI)= E12*CSA(I)*CSB(J)
            IF(JTYPE.GT.3) THEN
               D01P(JI)= E12*CSA(I)*CPB(J)
               IF(JTYPE.LE.5) THEN
                  IF(D01P(JI).NE.ZER) THEN
                     D00P(JI)= D00P(JI)/D01P(JI)
                  ENDIF
               ELSE
                  D10P(JI)= E12*CPA(I)*CSB(J)
                  D11P(JI)= E12*CPA(I)*CPB(J)
                  IF(D11P(JI).NE.ZER) THEN
                     TMP = ONE/D11P(JI)
                     D00P(JI)= D00P(JI)*TMP
                     D01P(JI)= D01P(JI)*TMP
                     D10P(JI)= D10P(JI)*TMP
                  ENDIF
               ENDIF
            ENDIF
  160       CONTINUE
  170 JI=JI+1
C
C BEGIN Q LOOP
C
      DO 190 K=1,NGC
         X03= EXC(K)
         DO 190 L=1,NGD
            X04= EXD(L)
            X34= X03+X04
            X43= ONE/X34
            Y03= X03*X43
            Y04= ONE-Y03
            Y34= Y03*X04
            R34Y34= R34*Y34
            IF(R34Y34.GT.CUX) GO TO 180
            E34= X43* EXP(-R34Y34)
            TST= E34*CMAXC(K)*CMAXD(L)
            IF(TST.LE.EI2) GO TO 180
            ISMLQ= 0
            IF(TST.LE.EI1) ISMLQ= 1
C
C CQX = COMPONENT OF CQ ALONG PENULTIMATE X-AXIS
C CQZ = COMPONENT OF CQ ALONG PENULTIMATE Z-AXIS
C
            CQ = RCD*Y04
            CQX= CQ*SING
            CQZ= CQ*COSG
C
C FIND COORDINATES OF Q RELATIVE TO AXES AT A
C QPR IS PERPENDICULAR FROM Q TO AB
C
            AQX= ACX+CQX
            AQX2= AQX*AQX
            AQXY= AQX*ACY
            AQZ= ACZ+CQZ
            QPS= AQX2+ACY2
C
            IF(JTYPE.NE.1) THEN
               SQ(0,0)= E34*CSC(K)*CSD(L)
               SQ(1,0)= E34*CSC(K)*CPD(L)
               SQ(0,1)= E34*CPC(K)*CSD(L)
               SQ(1,1)= E34*CPC(K)*CPD(L)
            ELSE
               SQ(0,0)= E34*CSC(K)*CSD(L)
            ENDIF
C
C USE SPECIAL FAST ROUTINE FOR INNER LOOPS FOR 0000 ... 1111
C
            CALL SP0S1S(JTYPE)
C
  180       CONTINUE
  190 CONTINUE
C
      IF(JTYPE.GT.1) THEN
C
CJMS  NOW, THE TRANSPOSE OF P TO BE USED FOR COMPUTATIONAL EFFICIENCY
C
         DO 195 J=1,2
            DO 195 I=J+1,3
               TMP= P(I,J)
               P(I,J)= P(J,I)
               P(J,I)= TMP
  195    CONTINUE
C
         QX= RCD*SING
         QZ= RCD*COSG
         CALL MCDAV(JTYPE,GPOPLE(1+NORG),QX,QZ)
         CALL R30S1S(JTYPE,GPOPLE(1+NORG),P)
      ELSE
         GPOPLE(1+NORG)= R00(1)
      ENDIF
  999 CONTINUE
      RETURN
      END
C
C*MODULE INT2B   *DECK MCDAV
      SUBROUTINE MCDAV(JTYPE,F,QX,QZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /KI3   / R00(25),R01(3,40),R02(6,26),R03(10,8),R04(15)
      DIMENSION F(0:3,0:3,0:3,0:3)
C
      IF(JTYPE.EQ.2)THEN
         F(0,0,0,0)= R00(1)
         F(1,0,0,0)= R01(1,1)+R00(2)*QX
         F(2,0,0,0)= R01(2,1)
         F(3,0,0,0)= R01(3,1)+R00(2)*QZ
      ELSEIF(JTYPE.EQ.3) THEN
         QXX= QX*QX
         QXZ= QX*QZ
         QZZ= QZ*QZ
         F(0,0,0,0)= R00(1)
         F(1,0,0,0)= R01(1,1)+R00(2)*QX
         F(2,0,0,0)= R01(2,1)
         F(3,0,0,0)= R01(3,1)+R00(2)*QZ
         F(0,1,0,0)= R01(1,2)+R00(3)*QX
         F(1,1,0,0)= R02(1,1)+R01(1,3)*QX+R01(1,4)*QX+R00(4)+R00(5)*QXX
         F(2,1,0,0)= R02(4,1)            +R01(2,4)*QX
         F(3,1,0,0)= R02(5,1)+R01(1,3)*QZ+R01(3,4)*QX       +R00(5)*QXZ
         F(0,2,0,0)= R01(2,2)
         F(1,2,0,0)= R02(4,1)+R01(2,3)*QX
         F(2,2,0,0)= R02(2,1)                        +R00(4)
         F(3,2,0,0)= R02(6,1)+R01(2,3)*QZ
         F(0,3,0,0)= R01(3,2)+R00(3)*QZ
         F(1,3,0,0)= R02(5,1)+R01(3,3)*QX+R01(1,4)*QZ       +R00(5)*QXZ
         F(2,3,0,0)= R02(6,1)            +R01(2,4)*QZ
         F(3,3,0,0)= R02(3,1)+R01(3,3)*QZ+R01(3,4)*QZ+R00(4)+R00(5)*QZZ
      ELSEIF(JTYPE.EQ.4) THEN
         F(0,0,0,0)= R00(1)
         F(1,0,0,0)= R01(1,1)+R00(2)*QX
         F(2,0,0,0)= R01(2,1)
         F(3,0,0,0)= R01(3,1)+R00(2)*QZ
         F(0,0,1,0)=-R01(1,2)
         F(1,0,1,0)=-R02(1,1)-R01(1,4)*QX
         F(2,0,1,0)=-R02(4,1)
         F(3,0,1,0)=-R02(5,1)-R01(1,4)*QZ
         F(0,0,2,0)=-R01(2,2)
         F(1,0,2,0)=-R02(4,1)-R01(2,4)*QX
         F(2,0,2,0)=-R02(2,1)
         F(3,0,2,0)=-R02(6,1)-R01(2,4)*QZ
         F(0,0,3,0)=-R01(3,2)+R00(3)
         F(1,0,3,0)=-R02(5,1)+R01(1,3)-R01(3,4)*QX+R00(4)*QX
         F(2,0,3,0)=-R02(6,1)+R01(2,3)
         F(3,0,3,0)=-R02(3,1)+R01(3,3)-R01(3,4)*QZ+R00(4)*QZ
      ELSEIF(JTYPE.EQ.5) THEN
         QXX= QX*QX
         QXZ= QX*QZ
         QZZ= QZ*QZ
         F(0,0,0,0)= R00(1)
         F(1,0,0,0)= R01(1,1)+R00(2)*QX
         F(2,0,0,0)= R01(2,1)
         F(3,0,0,0)= R01(3,1)+R00(2)*QZ
         F(0,1,0,0)= R01(1,2)+R00(3)*QX
         F(1,1,0,0)= R02(1,1)+R01(1,3)*QX+R01(1,4)*QX+R00(4)*QXX+R00(5)
         F(2,1,0,0)= R02(4,1)            +R01(2,4)*QX
         F(3,1,0,0)= R02(5,1)+R01(1,3)*QZ+R01(3,4)*QX+R00(4)*QXZ
         F(0,2,0,0)= R01(2,2)
         F(1,2,0,0)= R02(4,1)+R01(2,3)*QX
         F(2,2,0,0)= R02(2,1)                                   +R00(5)
         F(3,2,0,0)= R02(6,1)+R01(2,3)*QZ
         F(0,3,0,0)= R01(3,2)+R00(3)*QZ
         F(1,3,0,0)= R02(5,1)+R01(3,3)*QX+R01(1,4)*QZ+R00(4)*QXZ
         F(2,3,0,0)= R02(6,1)            +R01(2,4)*QZ
         F(3,3,0,0)= R02(3,1)+R01(3,3)*QZ+R01(3,4)*QZ+R00(4)*QZZ+R00(5)
         F(0,0,1,0)=-R01(1,9)
         F(1,0,1,0)=-R02(1,3)-R01(1,10)*QX
         F(2,0,1,0)=-R02(4,3)
         F(3,0,1,0)=-R02(5,3)-R01(1,10)*QZ
         F(0,1,1,0)=-R02(1,4)-R01(1,11)*QX
         F(1,1,1,0)=-R03(1,1)-R02(1,5)*QX-R02(1,6)*QX-R01(1,12)*QXX
     *                                               -R01(1,13)
         F(2,1,1,0)=-R03(2,1)            -R02(4,6)*QX
         F(3,1,1,0)=-R03(3,1)-R02(1,5)*QZ-R02(5,6)*QX-R01(1,12)*QXZ
         F(0,2,1,0)=-R02(4,4)
         F(1,2,1,0)=-R03(2,1)-R02(4,5)*QX
         F(2,2,1,0)=-R03(4,1)                        -R01(1,13)
         F(3,2,1,0)=-R03(5,1)-R02(4,5)*QZ
         F(0,3,1,0)=-R02(5,4)-R01(1,11)*QZ
         F(1,3,1,0)=-R03(3,1)-R02(5,5)*QX-R02(1,6)*QZ-R01(1,12)*QXZ
         F(2,3,1,0)=-R03(5,1)            -R02(4,6)*QZ
         F(3,3,1,0)=-R03(6,1)-R02(5,5)*QZ-R02(5,6)*QZ-R01(1,12)*QZZ
     *                                               -R01(1,13)
         F(0,0,2,0)=-R01(2,9)
         F(1,0,2,0)=-R02(4,3)-R01(2,10)*QX
         F(2,0,2,0)=-R02(2,3)
         F(3,0,2,0)=-R02(6,3)-R01(2,10)*QZ
         F(0,1,2,0)=-R02(4,4)-R01(2,11)*QX
         F(1,1,2,0)=-R03(2,1)-R02(4,5)*QX-R02(4,6)*QX-R01(2,12)*QXX
     *                                               -R01(2,13)
         F(2,1,2,0)=-R03(4,1)            -R02(2,6)*QX
         F(3,1,2,0)=-R03(5,1)-R02(4,5)*QZ-R02(6,6)*QX-R01(2,12)*QXZ
         F(0,2,2,0)=-R02(2,4)
         F(1,2,2,0)=-R03(4,1)-R02(2,5)*QX
         F(2,2,2,0)=-R03(7,1)                        -R01(2,13)
         F(3,2,2,0)=-R03(8,1)-R02(2,5)*QZ
         F(0,3,2,0)=-R02(6,4)-R01(2,11)*QZ
         F(1,3,2,0)=-R03(5,1)-R02(6,5)*QX-R02(4,6)*QZ-R01(2,12)*QXZ
         F(2,3,2,0)=-R03(8,1)            -R02(2,6)*QZ
         F(3,3,2,0)=-R03(9,1)-R02(6,5)*QZ-R02(6,6)*QZ-R01(2,12)*QZZ
     *                                               -R01(2,13)
         F(0,0,3,0)=-R01(3,9)+R00(6)
         F(1,0,3,0)=-R02(5,3)+R01(1,5)-R01(3,10)*QX+R00(7)*QX
         F(2,0,3,0)=-R02(6,3)+R01(2,5)
         F(3,0,3,0)=-R02(3,3)+R01(3,5)-R01(3,10)*QZ+R00(7)*QZ
         F(0,1,3,0)=-R02(5,4)+R01(1,6)-R01(3,11)*QX+R00(8)*QX
         F(1,1,3,0)=-R03(3,1)+R02(1,2)-R02(5,5)*QX-R02(5,6)*QX
     *              +R01(1,7)*QX+R01(1,8)*QX-R01(3,12)*QXX-R01(3,13)
     *              +R00(9)*QXX+R00(10)
         F(2,1,3,0)=-R03(5,1)+R02(4,2)-R02(6,6)*QX+R01(2,8)*QX
         F(3,1,3,0)=-R03(6,1)+R02(5,2)-R02(5,5)*QZ-R02(3,6)*QX
     *              +R01(1,7)*QZ+R01(3,8)*QX-R01(3,12)*QXZ+R00(9)*QXZ
         F(0,2,3,0)=-R02(6,4)+R01(2,6)
         F(1,2,3,0)=-R03(5,1)+R02(4,2)-R02(6,5)*QX+R01(2,7)*QX
         F(2,2,3,0)=-R03(8,1)+R02(2,2)-R01(3,13)+R00(10)
         F(3,2,3,0)=-R03(9,1)+R02(6,2)-R02(6,5)*QZ+R01(2,7)*QZ
         F(0,3,3,0)=-R02(3,4)+R01(3,6)-R01(3,11)*QZ+R00(8)*QZ
         F(1,3,3,0)=-R03(6,1)+R02(5,2)-R02(3,5)*QX-R02(5,6)*QZ
     *              +R01(3,7)*QX+R01(1,8)*QZ-R01(3,12)*QXZ+R00(9)*QXZ
         F(2,3,3,0)=-R03(9,1)+R02(6,2)-R02(6,6)*QZ+R01(2,8)*QZ
         F(3,3,3,0)=-R03(10,1)+R02(3,2)-R02(3,5)*QZ-R02(3,6)*QZ
     *              +R01(3,7)*QZ+R01(3,8)*QZ-R01(3,12)*QZZ-R01(3,13)
     *              +R00(9)*QZZ+R00(10)
      ELSEIF(JTYPE.EQ.6) THEN
         QXX= QX*QX
         QXZ= QX*QZ
         QZZ= QZ*QZ
         F(0,0,0,0)= R00(1)
         F(1,0,0,0)= R01(1,1)+R00(2)*QX
         F(2,0,0,0)= R01(2,1)
         F(3,0,0,0)= R01(3,1)+R00(2)*QZ
         F(0,1,0,0)= R01(1,2)+R00(3)*QX
         F(1,1,0,0)= R02(1,1)+R01(1,3)*QX+R01(1,4)*QX+R00(4)*QXX+R00(5)
         F(2,1,0,0)= R02(4,1)            +R01(2,4)*QX
         F(3,1,0,0)= R02(5,1)+R01(1,3)*QZ+R01(3,4)*QX+R00(4)*QXZ
         F(0,2,0,0)= R01(2,2)
         F(1,2,0,0)= R02(4,1)+R01(2,3)*QX
         F(2,2,0,0)= R02(2,1)                                   +R00(5)
         F(3,2,0,0)= R02(6,1)+R01(2,3)*QZ
         F(0,3,0,0)= R01(3,2)+R00(3)*QZ
         F(1,3,0,0)= R02(5,1)+R01(3,3)*QX+R01(1,4)*QZ+R00(4)*QXZ
         F(2,3,0,0)= R02(6,1)            +R01(2,4)*QZ
         F(3,3,0,0)= R02(3,1)+R01(3,3)*QZ+R01(3,4)*QZ+R00(4)*QZZ+R00(5)
         F(0,0,1,0)=-R01(1,21)
         F(1,0,1,0)=-R02(1,6)-R01(1,22)*QX
         F(2,0,1,0)=-R02(4,6)
         F(3,0,1,0)=-R02(5,6)-R01(1,22)*QZ
         F(0,1,1,0)=-R02(1,7)-R01(1,23)*QX
         F(1,1,1,0)=-R03(1,2)-R02(1,8)*QX-R02(1,9)*QX-R01(1,24)*QXX
     *                                               -R01(1,25)
         F(2,1,1,0)=-R03(2,2)            -R02(4,9)*QX
         F(3,1,1,0)=-R03(3,2)-R02(1,8)*QZ-R02(5,9)*QX-R01(1,24)*QXZ
         F(0,2,1,0)=-R02(4,7)
         F(1,2,1,0)=-R03(2,2)-R02(4,8)*QX
         F(2,2,1,0)=-R03(4,2)                        -R01(1,25)
         F(3,2,1,0)=-R03(5,2)-R02(4,8)*QZ
         F(0,3,1,0)=-R02(5,7)-R01(1,23)*QZ
         F(1,3,1,0)=-R03(3,2)-R02(5,8)*QX-R02(1,9)*QZ-R01(1,24)*QXZ
         F(2,3,1,0)=-R03(5,2)            -R02(4,9)*QZ
         F(3,3,1,0)=-R03(6,2)-R02(5,8)*QZ-R02(5,9)*QZ-R01(1,24)*QZZ
     *                                               -R01(1,25)
C
         F(0,0,2,0)=-R01(2,21)
         F(1,0,2,0)=-R02(4,6)-R01(2,22)*QX
         F(2,0,2,0)=-R02(2,6)
         F(3,0,2,0)=-R02(6,6)-R01(2,22)*QZ
         F(0,1,2,0)=-R02(4,7)-R01(2,23)*QX
         F(1,1,2,0)=-R03(2,2)-R02(4,8)*QX-R02(4,9)*QX-R01(2,24)*QXX
     *                                               -R01(2,25)
         F(2,1,2,0)=-R03(4,2)            -R02(2,9)*QX
         F(3,1,2,0)=-R03(5,2)-R02(4,8)*QZ-R02(6,9)*QX-R01(2,24)*QXZ
         F(0,2,2,0)=-R02(2,7)
         F(1,2,2,0)=-R03(4,2)-R02(2,8)*QX
         F(2,2,2,0)=-R03(7,2)                        -R01(2,25)
         F(3,2,2,0)=-R03(8,2)-R02(2,8)*QZ
         F(0,3,2,0)=-R02(6,7)-R01(2,23)*QZ
         F(1,3,2,0)=-R03(5,2)-R02(6,8)*QX-R02(4,9)*QZ-R01(2,24)*QXZ
         F(2,3,2,0)=-R03(8,2)            -R02(2,9)*QZ
         F(3,3,2,0)=-R03(9,2)-R02(6,8)*QZ-R02(6,9)*QZ-R01(2,24)*QZZ
     *                                               -R01(2,25)
         F(0,0,3,0)=-R01(3,21)+R00(6)
         F(1,0,3,0)=-R02(5,6)+R01(1,5)-R01(3,22)*QX+R00(7)*QX
         F(2,0,3,0)=-R02(6,6)+R01(2,5)
         F(3,0,3,0)=-R02(3,6)+R01(3,5)-R01(3,22)*QZ+R00(7)*QZ
         F(0,1,3,0)=-R02(5,7)+R01(1,6)-R01(3,23)*QX+R00(8)*QX
         F(1,1,3,0)=-R03(3,2)+R02(1,2)-R02(5,8)*QX-R02(5,9)*QX
     *              +R01(1,7)*QX+R01(1,8)*QX-R01(3,24)*QXX-R01(3,25)
     *              +R00(9)*QXX+R00(10)
         F(2,1,3,0)=-R03(5,2)+R02(4,2)-R02(6,9)*QX+R01(2,8)*QX
         F(3,1,3,0)=-R03(6,2)+R02(5,2)-R02(5,8)*QZ-R02(3,9)*QX
     *              +R01(1,7)*QZ+R01(3,8)*QX-R01(3,24)*QXZ+R00(9)*QXZ
         F(0,2,3,0)=-R02(6,7)+R01(2,6)
         F(1,2,3,0)=-R03(5,2)+R02(4,2)-R02(6,8)*QX+R01(2,7)*QX
         F(2,2,3,0)=-R03(8,2)+R02(2,2)-R01(3,25)+R00(10)
         F(3,2,3,0)=-R03(9,2)+R02(6,2)-R02(6,8)*QZ+R01(2,7)*QZ
         F(0,3,3,0)=-R02(3,7)+R01(3,6)-R01(3,23)*QZ+R00(8)*QZ
         F(1,3,3,0)=-R03(6,2)+R02(5,2)-R02(3,8)*QX-R02(5,9)*QZ
     *              +R01(3,7)*QX+R01(1,8)*QZ-R01(3,24)*QXZ+R00(9)*QXZ
         F(2,3,3,0)=-R03(9,2)+R02(6,2)-R02(6,9)*QZ+R01(2,8)*QZ
         F(3,3,3,0)=-R03(10,2)+R02(3,2)-R02(3,8)*QZ-R02(3,9)*QZ
     *              +R01(3,7)*QZ+R01(3,8)*QZ-R01(3,24)*QZZ-R01(3,25)
     *              +R00(9)*QZZ+R00(10)
C
         F(0,0,0,1)=-R01(1,26)
         F(1,0,0,1)=-R02(1,10)-R01(1,27)*QX
         F(2,0,0,1)=-R02(4,10)
         F(3,0,0,1)=-R02(5,10)-R01(1,27)*QZ
         F(0,1,0,1)=-R02(1,11)-R01(1,28)*QX
         F(1,1,0,1)=-R03(1,3)-R02(1,12)*QX-R02(1,13)*QX-R01(1,29)*QXX
     *                                                 -R01(1,30)
         F(2,1,0,1)=-R03(2,3)             -R02(4,13)*QX
         F(3,1,0,1)=-R03(3,3)-R02(1,12)*QZ-R02(5,13)*QX-R01(1,29)*QXZ
         F(0,2,0,1)=-R02(4,11)
         F(1,2,0,1)=-R03(2,3)-R02(4,12)*QX
         F(2,2,0,1)=-R03(4,3)                          -R01(1,30)
         F(3,2,0,1)=-R03(5,3)-R02(4,12)*QZ
         F(0,3,0,1)=-R02(5,11)-R01(1,28)*QZ
         F(1,3,0,1)=-R03(3,3)-R02(5,12)*QX-R02(1,13)*QZ-R01(1,29)*QXZ
         F(2,3,0,1)=-R03(5,3)             -R02(4,13)*QZ
         F(3,3,0,1)=-R03(6,3)-R02(5,12)*QZ-R02(5,13)*QZ-R01(1,29)*QZZ
     *                                                 -R01(1,30)
         F(0,0,1,1)= R02(1,22)+R00(21)
         F(1,0,1,1)= R03(1,5)+R02(1,23)*QX+R01(1,17)+R00(22)*QX
         F(2,0,1,1)= R03(2,5)             +R01(2,17)
         F(3,0,1,1)= R03(3,5)+R02(1,23)*QZ+R01(3,17)+R00(22)*QZ
         F(0,1,1,1)= R03(1,6)+R02(1,24)*QX+R01(1,18)+R00(23)*QX
         F(1,1,1,1)= R04( 1)+R03(1,7)*QX+R03(1,8)*QX+R02(1,5)
     *              +R02(1,25)*QXX+R02(1,26)+R01(1,19)*QX+R01(1,20)*QX
     *              +R00(24)*QXX+R00(25)
         F(2,1,1,1)= R04( 2)+R03(2,8)*QX+R02(4,5)+R01(2,20)*QX
         F(3,1,1,1)= R04( 3)+R03(1,7)*QZ+R03(3,8)*QX+R02(5,5)
     *              +R02(1,25)*QXZ+R01(1,19)*QZ+R01(3,20)*QX+R00(24)*QXZ
         F(0,2,1,1)= R03(2,6)             +R01(2,18)
         F(1,2,1,1)= R04( 2)+R03(2,7)*QX+R02(4,5)+R01(2,19)*QX
         F(2,2,1,1)= R04( 4)+R02(2,5)+R02(1,26)+R00(25)
         F(3,2,1,1)= R04( 5)+R03(2,7)*QZ+R02(6,5)+R01(2,19)*QZ
         F(0,3,1,1)= R03(3,6)+R02(1,24)*QZ+R01(3,18)+R00(23)*QZ
         F(1,3,1,1)= R04( 3)+R03(3,7)*QX+R03(1,8)*QZ+R02(5,5)
     *              +R02(1,25)*QXZ+R01(3,19)*QX+R01(1,20)*QZ+R00(24)*QXZ
         F(2,3,1,1)= R04( 5)+R03(2,8)*QZ+R02(6,5)+R01(2,20)*QZ
         F(3,3,1,1)= R04( 6)+R03(3,7)*QZ+R03(3,8)*QZ+R02(3,5)
     *              +R02(1,25)*QZZ+R02(1,26)+R01(3,19)*QZ+R01(3,20)*QZ
     *              +R00(24)*QZZ+R00(25)
         F(0,0,2,1)= R02(4,22)
         F(1,0,2,1)= R03(2,5)+R02(4,23)*QX
         F(2,0,2,1)= R03(4,5)
         F(3,0,2,1)= R03(5,5)+R02(4,23)*QZ
         F(0,1,2,1)= R03(2,6)+R02(4,24)*QX
         F(1,1,2,1)= R04( 2)+R03(2,7)*QX+R03(2,8)*QX+R02(4,25)*QXX
     *                                              +R02(4,26)
         F(2,1,2,1)= R04( 4)            +R03(4,8)*QX
         F(3,1,2,1)= R04( 5)+R03(2,7)*QZ+R03(5,8)*QX+R02(4,25)*QXZ
         F(0,2,2,1)= R03(4,6)
         F(1,2,2,1)= R04( 4)+R03(4,7)*QX
         F(2,2,2,1)= R04( 7)                        +R02(4,26)
         F(3,2,2,1)= R04( 8)+R03(4,7)*QZ
         F(0,3,2,1)= R03(5,6)+R02(4,24)*QZ
         F(1,3,2,1)= R04( 5)+R03(5,7)*QX+R03(2,8)*QZ+R02(4,25)*QXZ
         F(2,3,2,1)= R04( 8)            +R03(4,8)*QZ
         F(3,3,2,1)= R04( 9)+R03(5,7)*QZ+R03(5,8)*QZ+R02(4,25)*QZZ
     *                                              +R02(4,26)
         F(0,0,3,1)= R02(5,22)-R01(1,31)
         F(1,0,3,1)= R03(3,5)+R02(5,23)*QX-R02(1,14)-R01(1,32)*QX
         F(2,0,3,1)= R03(5,5)             -R02(4,14)
         F(3,0,3,1)= R03(6,5)+R02(5,23)*QZ-R02(5,14)-R01(1,32)*QZ
         F(0,1,3,1)= R03(3,6)+R02(5,24)*QX-R02(1,15)-R01(1,33)*QX
         F(1,1,3,1)= R04( 3)-R03(1,4)+R03(3,7)*QX+R03(3,8)*QX
     *              -R02(1,16)*QX-R02(1,17)*QX+R02(5,25)*QXX+R02(5,26)
     *              -R01(1,34)*QXX-R01(1,35)
         F(2,1,3,1)= R04( 5)-R03(2,4)+R03(5,8)*QX-R02(4,17)*QX
         F(3,1,3,1)= R04( 6)-R03(3,4)+R03(3,7)*QZ+R03(6,8)*QX
     *              -R02(1,16)*QZ-R02(5,17)*QX+R02(5,25)*QXZ
     *              -R01(1,34)*QXZ
         F(0,2,3,1)= R03(5,6)             -R02(4,15)
         F(1,2,3,1)= R04( 5)-R03(2,4)+R03(5,7)*QX-R02(4,16)*QX
         F(2,2,3,1)= R04( 8)-R03(4,4)+R02(5,26)-R01(1,35)
         F(3,2,3,1)= R04( 9)-R03(5,4)+R03(5,7)*QZ-R02(4,16)*QZ
         F(0,3,3,1)= R03(6,6)+R02(5,24)*QZ-R02(5,15)-R01(1,33)*QZ
         F(1,3,3,1)= R04( 6)-R03(3,4)+R03(6,7)*QX+R03(3,8)*QZ
     *              -R02(5,16)*QX-R02(1,17)*QZ+R02(5,25)*QXZ
     *              -R01(1,34)*QXZ
         F(2,3,3,1)= R04( 9)-R03(5,4)+R03(5,8)*QZ-R02(4,17)*QZ
         F(3,3,3,1)= R04(10)-R03(6,4)+R03(6,7)*QZ+R03(6,8)*QZ
     *              -R02(5,16)*QZ-R02(5,17)*QZ+R02(5,25)*QZZ+R02(5,26)
     *              -R01(1,34)*QZZ-R01(1,35)
C
         F(0,0,0,2)=-R01(2,26)
         F(1,0,0,2)=-R02(4,10)-R01(2,27)*QX
         F(2,0,0,2)=-R02(2,10)
         F(3,0,0,2)=-R02(6,10)-R01(2,27)*QZ
         F(0,1,0,2)=-R02(4,11)-R01(2,28)*QX
         F(1,1,0,2)=-R03(2,3)-R02(4,12)*QX-R02(4,13)*QX-R01(2,29)*QXX
     *                                                 -R01(2,30)
         F(2,1,0,2)=-R03(4,3)             -R02(2,13)*QX
         F(3,1,0,2)=-R03(5,3)-R02(4,12)*QZ-R02(6,13)*QX-R01(2,29)*QXZ
         F(0,2,0,2)=-R02(2,11)
         F(1,2,0,2)=-R03(4,3)-R02(2,12)*QX
         F(2,2,0,2)=-R03(7,3)                          -R01(2,30)
         F(3,2,0,2)=-R03(8,3)-R02(2,12)*QZ
         F(0,3,0,2)=-R02(6,11)-R01(2,28)*QZ
         F(1,3,0,2)=-R03(5,3)-R02(6,12)*QX-R02(4,13)*QZ-R01(2,29)*QXZ
         F(2,3,0,2)=-R03(8,3)             -R02(2,13)*QZ
         F(3,3,0,2)=-R03(9,3)-R02(6,12)*QZ-R02(6,13)*QZ-R01(2,29)*QZZ
     *                                                 -R01(2,30)
         F(0,0,1,2)= R02(4,22)
         F(1,0,1,2)= R03(2,5)+R02(4,23)*QX
         F(2,0,1,2)= R03(4,5)
         F(3,0,1,2)= R03(5,5)+R02(4,23)*QZ
         F(0,1,1,2)= R03(2,6)+R02(4,24)*QX
         F(1,1,1,2)= R04( 2)+R03(2,7)*QX+R03(2,8)*QX+R02(4,25)*QXX
     *                                              +R02(4,26)
         F(2,1,1,2)= R04( 4)            +R03(4,8)*QX
         F(3,1,1,2)= R04( 5)+R03(2,7)*QZ+R03(5,8)*QX+R02(4,25)*QXZ
         F(0,2,1,2)= R03(4,6)
         F(1,2,1,2)= R04( 4)+R03(4,7)*QX
         F(2,2,1,2)= R04( 7)                        +R02(4,26)
         F(3,2,1,2)= R04( 8)+R03(4,7)*QZ
         F(0,3,1,2)= R03(5,6)+R02(4,24)*QZ
         F(1,3,1,2)= R04( 5)+R03(5,7)*QX+R03(2,8)*QZ+R02(4,25)*QXZ
         F(2,3,1,2)= R04( 8)            +R03(4,8)*QZ
         F(3,3,1,2)= R04( 9)+R03(5,7)*QZ+R03(5,8)*QZ+R02(4,25)*QZZ
     *                                              +R02(4,26)
         F(0,0,2,2)= R02(2,22)+R00(21)
         F(1,0,2,2)= R03(4,5)+R02(2,23)*QX+R01(1,17)+R00(22)*QX
         F(2,0,2,2)= R03(7,5)             +R01(2,17)
         F(3,0,2,2)= R03(8,5)+R02(2,23)*QZ+R01(3,17)+R00(22)*QZ
         F(0,1,2,2)= R03(4,6)+R02(2,24)*QX+R01(1,18)+R00(23)*QX
         F(1,1,2,2)= R04( 4)+R03(4,7)*QX+R03(4,8)*QX+R02(1,5)
     *              +R02(2,25)*QXX+R02(2,26)+R01(1,19)*QX+R01(1,20)*QX
     *              +R00(24)*QXX+R00(25)
         F(2,1,2,2)= R04( 7)+R03(7,8)*QX+R02(4,5)+R01(2,20)*QX
         F(3,1,2,2)= R04( 8)+R03(4,7)*QZ+R03(8,8)*QX+R02(5,5)
     *              +R02(2,25)*QXZ+R01(1,19)*QZ+R01(3,20)*QX+R00(24)*QXZ
         F(0,2,2,2)= R03(7,6)             +R01(2,18)
         F(1,2,2,2)= R04( 7)+R03(7,7)*QX+R02(4,5)+R01(2,19)*QX
         F(2,2,2,2)= R04(11)+R02(2,5)+R02(2,26)+R00(25)
         F(3,2,2,2)= R04(12)+R03(7,7)*QZ+R02(6,5)+R01(2,19)*QZ
         F(0,3,2,2)= R03(8,6)+R02(2,24)*QZ+R01(3,18)+R00(23)*QZ
         F(1,3,2,2)= R04( 8)+R03(8,7)*QX+R03(4,8)*QZ+R02(5,5)
     *              +R02(2,25)*QXZ+R01(3,19)*QX+R01(1,20)*QZ+R00(24)*QXZ
         F(2,3,2,2)= R04(12)+R03(7,8)*QZ+R02(6,5)+R01(2,20)*QZ
         F(3,3,2,2)= R04(13)+R03(8,7)*QZ+R03(8,8)*QZ+R02(3,5)
     *              +R02(2,25)*QZZ+R02(2,26)+R01(3,19)*QZ+R01(3,20)*QZ
     *              +R00(24)*QZZ+R00(25)
         F(0,0,3,2)= R02(6,22)-R01(2,31)
         F(1,0,3,2)= R03(5,5)-R02(4,14)+R02(6,23)*QX-R01(2,32)*QX
         F(2,0,3,2)= R03(8,5)-R02(2,14)
         F(3,0,3,2)= R03(9,5)-R02(6,14)+R02(6,23)*QZ-R01(2,32)*QZ
         F(0,1,3,2)= R03(5,6)-R02(4,15)+R02(6,24)*QX-R01(2,33)*QX
         F(1,1,3,2)= R04( 5)+R03(5,7)*QX+R03(5,8)*QX-R03(2,4)
     *              -R02(4,16)*QX-R02(4,17)*QX+R02(6,25)*QXX+R02(6,26)
     *              -R01(2,34)*QXX-R01(2,35)
         F(2,1,3,2)= R04( 8)-R03(4,4)+R03(8,8)*QX-R02(2,17)*QX
         F(3,1,3,2)= R04( 9)-R03(5,4)+R03(5,7)*QZ+R03(9,8)*QX
     *              -R02(4,16)*QZ-R02(6,17)*QX+R02(6,25)*QXZ
     *              -R01(2,34)*QXZ
         F(0,2,3,2)= R03(8,6)-R02(2,15)
         F(1,2,3,2)= R04( 8)-R03(4,4)+R03(8,7)*QX-R02(2,16)*QX
         F(2,2,3,2)= R04(12)-R03(7,4)+R02(6,26)-R01(2,35)
         F(3,2,3,2)= R04(13)-R03(8,4)+R03(8,7)*QZ-R02(2,16)*QZ
         F(0,3,3,2)= R03(9,6)-R02(6,15)+R02(6,24)*QZ-R01(2,33)*QZ
         F(1,3,3,2)= R04( 9)-R03(5,4)+R03(9,7)*QX+R03(5,8)*QZ
     *              -R02(6,16)*QX-R02(4,17)*QZ+R02(6,25)*QXZ
     *              -R01(2,34)*QXZ
         F(2,3,3,2)= R04(13)-R03(8,4)+R03(8,8)*QZ-R02(2,17)*QZ
         F(3,3,3,2)= R04(14)-R03(9,4)+R03(9,7)*QZ+R03(9,8)*QZ
     *              -R02(6,16)*QZ-R02(6,17)*QZ+R02(6,25)*QZZ+R02(6,26)
     *              -R01(2,34)*QZZ-R01(2,35)
C
         F(0,0,0,3)=-R01(3,26)+R00(11)
         F(1,0,0,3)=-R02(5,10)+R01(1,9)-R01(3,27)*QX+R00(12)*QX
         F(2,0,0,3)=-R02(6,10)+R01(2,9)
         F(3,0,0,3)=-R02(3,10)+R01(3,9)-R01(3,27)*QZ+R00(12)*QZ
         F(0,1,0,3)=-R02(5,11)+R01(1,10)-R01(3,28)*QX+R00(13)*QX
         F(1,1,0,3)=-R03(3,3)+R02(1,3)-R02(5,12)*QX-R02(5,13)*QX
     *              +R01(1,11)*QX+R01(1,12)*QX-R01(3,29)*QXX-R01(3,30)
     *              +R00(14)*QXX+R00(15)
         F(2,1,0,3)=-R03(5,3)+R02(4,3)-R02(6,13)*QX+R01(2,12)*QX
         F(3,1,0,3)=-R03(6,3)+R02(5,3)-R02(5,12)*QZ-R02(3,13)*QX
     *              +R01(1,11)*QZ+R01(3,12)*QX-R01(3,29)*QXZ+R00(14)*QXZ
         F(0,2,0,3)=-R02(6,11)+R01(2,10)
         F(1,2,0,3)=-R03(5,3)+R02(4,3)-R02(6,12)*QX+R01(2,11)*QX
         F(2,2,0,3)=-R03(8,3)+R02(2,3)-R01(3,30)+R00(15)
         F(3,2,0,3)=-R03(9,3)+R02(6,3)-R02(6,12)*QZ+R01(2,11)*QZ
         F(0,3,0,3)=-R02(3,11)+R01(3,10)-R01(3,28)*QZ+R00(13)*QZ
         F(1,3,0,3)=-R03(6,3)+R02(5,3)-R02(3,12)*QX-R02(5,13)*QZ
     *              +R01(3,11)*QX+R01(1,12)*QZ-R01(3,29)*QXZ+R00(14)*QXZ
         F(2,3,0,3)=-R03(9,3)+R02(6,3)-R02(6,13)*QZ+R01(2,12)*QZ
         F(3,3,0,3)=-R03(10,3)+R02(3,3)-R02(3,12)*QZ-R02(3,13)*QZ
     *              +R01(3,11)*QZ+R01(3,12)*QZ-R01(3,29)*QZZ-R01(3,30)
     *              +R00(14)*QZZ+R00(15)
         F(0,0,1,3)= R02(5,22)-R01(1,36)
         F(1,0,1,3)= R03(3,5)-R02(1,18)+R02(5,23)*QX-R01(1,37)*QX
         F(2,0,1,3)= R03(5,5)-R02(4,18)
         F(3,0,1,3)= R03(6,5)-R02(5,18)+R02(5,23)*QZ-R01(1,37)*QZ
         F(0,1,1,3)= R03(3,6)-R02(1,19)+R02(5,24)*QX-R01(1,38)*QX
         F(1,1,1,3)= R04( 3)-R03(1,1)+R03(3,7)*QX+R03(3,8)*QX
     *              -R02(1,20)*QX-R02(1,21)*QX+R02(5,25)*QXX+R02(5,26)
     *              -R01(1,39)*QXX-R01(1,40)
         F(2,1,1,3)= R04( 5)-R03(2,1)+R03(5,8)*QX-R02(4,21)*QX
         F(3,1,1,3)= R04( 6)-R03(3,1)+R03(3,7)*QZ+R03(6,8)*QX
     *              -R02(1,20)*QZ-R02(5,21)*QX+R02(5,25)*QXZ
     *              -R01(1,39)*QXZ
         F(0,2,1,3)= R03(5,6)-R02(4,19)
         F(1,2,1,3)= R04( 5)-R03(2,1)+R03(5,7)*QX-R02(4,20)*QX
         F(2,2,1,3)= R04( 8)-R03(4,1)+R02(5,26)-R01(1,40)
         F(3,2,1,3)= R04( 9)-R03(5,1)+R03(5,7)*QZ-R02(4,20)*QZ
         F(0,3,1,3)= R03(6,6)-R02(5,19)+R02(5,24)*QZ-R01(1,38)*QZ
         F(1,3,1,3)= R04( 6)-R03(3,1)+R03(6,7)*QX+R03(3,8)*QZ
     *              -R02(5,20)*QX-R02(1,21)*QZ+R02(5,25)*QXZ
     *              -R01(1,39)*QXZ
         F(2,3,1,3)= R04( 9)-R03(5,1)+R03(5,8)*QZ-R02(4,21)*QZ
         F(3,3,1,3)= R04(10)-R03(6,1)+R03(6,7)*QZ+R03(6,8)*QZ
     *              -R02(5,20)*QZ-R02(5,21)*QZ+R02(5,25)*QZZ+R02(5,26)
     *              -R01(1,39)*QZZ-R01(1,40)
         F(0,0,2,3)= R02(6,22)-R01(2,36)
         F(1,0,2,3)= R03(5,5)-R02(4,18)+R02(6,23)*QX-R01(2,37)*QX
         F(2,0,2,3)= R03(8,5)-R02(2,18)
         F(3,0,2,3)= R03(9,5)-R02(6,18)+R02(6,23)*QZ-R01(2,37)*QZ
         F(0,1,2,3)= R03(5,6)-R02(4,19)+R02(6,24)*QX-R01(2,38)*QX
         F(1,1,2,3)= R04( 5)-R03(2,1)+R03(5,7)*QX+R03(5,8)*QX
     *              -R02(4,20)*QX-R02(4,21)*QX+R02(6,25)*QXX+R02(6,26)
     *              -R01(2,39)*QXX-R01(2,40)
         F(2,1,2,3)= R04( 8)-R03(4,1)+R03(8,8)*QX-R02(2,21)*QX
         F(3,1,2,3)= R04( 9)-R03(5,1)+R03(5,7)*QZ+R03(9,8)*QX
     *              -R02(4,20)*QZ-R02(6,21)*QX+R02(6,25)*QXZ
     *              -R01(2,39)*QXZ
         F(0,2,2,3)= R03(8,6)-R02(2,19)
         F(1,2,2,3)= R04( 8)-R03(4,1)+R03(8,7)*QX-R02(2,20)*QX
         F(2,2,2,3)= R04(12)-R03(7,1)+R02(6,26)-R01(2,40)
         F(3,2,2,3)= R04(13)-R03(8,1)+R03(8,7)*QZ-R02(2,20)*QZ
         F(0,3,2,3)= R03(9,6)-R02(6,19)+R02(6,24)*QZ-R01(2,38)*QZ
         F(1,3,2,3)= R04( 9)-R03(5,1)+R03(9,7)*QX+R03(5,8)*QZ
     *              -R02(6,20)*QX-R02(4,21)*QZ+R02(6,25)*QXZ
     *              -R01(2,39)*QXZ
         F(2,3,2,3)= R04(13)-R03(8,1)+R03(8,8)*QZ-R02(2,21)*QZ
         F(3,3,2,3)= R04(14)-R03(9,1)+R03(9,7)*QZ+R03(9,8)*QZ
     *              -R02(6,20)*QZ-R02(6,21)*QZ+R02(6,25)*QZZ+R02(6,26)
     *              -R01(2,39)*QZZ-R01(2,40)
         F(0,0,3,3)= R02(3,22)-R01(3,31)-R01(3,36)+R00(16)+R00(21)
         F(1,0,3,3)= R03(6,5)-R02(5,14)-R02(5,18)+R02(3,23)*QX
     *              +R01(1,13)+R01(1,17)-R01(3,32)*QX-R01(3,37)*QX
     *              +R00(17)*QX+R00(22)*QX
         F(2,0,3,3)= R03(9,5)-R02(6,14)-R02(6,18)+R01(2,13)+R01(2,17)
         F(3,0,3,3)= R03(10,5)-R02(3,14)-R02(3,18)+R02(3,23)*QZ
     *              +R01(3,13)+R01(3,17)-R01(3,32)*QZ-R01(3,37)*QZ
     *              +R00(17)*QZ+R00(22)*QZ
         F(0,1,3,3)= R03(6,6)-R02(5,15)-R02(5,19)+R02(3,24)*QX
     *              +R01(1,14)+R01(1,18)-R01(3,33)*QX-R01(3,38)*QX
     *              +R00(18)*QX+R00(23)*QX
         F(1,1,3,3)= R04( 6)-R03(3,1)-R03(3,4)+R03(6,7)*QX+R03(6,8)*QX
     *              +R02(1,4)+R02(1,5)-R02(5,16)*QX-R02(5,17)*QX
     *              -R02(5,20)*QX-R02(5,21)*QX+R02(3,25)*QXX+R02(3,26)
     *              +R01(1,15)*QX+R01(1,16)*QX+R01(1,19)*QX+R01(1,20)*QX
     *              -R01(3,34)*QXX-R01(3,35)-R01(3,39)*QXX-R01(3,40)
     *              +R00(19)*QXX+R00(20)+R00(24)*QXX+R00(25)
         F(2,1,3,3)= R04( 9)-R03(5,1)-R03(5,4)+R03(9,8)*QX
     *              +R02(4,4)+R02(4,5)-R02(6,17)*QX-R02(6,21)*QX
     *              +R01(2,16)*QX+R01(2,20)*QX
         F(3,1,3,3)= R04(10)-R03(6,1)-R03(6,4)+R03(6,7)*QZ+R03(10,8)*QX
     *              +R02(5,4)+R02(5,5)-R02(5,16)*QZ-R02(3,17)*QX
     *              -R02(5,20)*QZ-R02(3,21)*QX+R02(3,25)*QXZ
     *              +R01(1,15)*QZ+R01(3,16)*QX+R01(1,19)*QZ+R01(3,20)*QX
     *              -R01(3,34)*QXZ-R01(3,39)*QXZ+R00(19)*QXZ+R00(24)*QXZ
         F(0,2,3,3)= R03(9,6)-R02(6,15)-R02(6,19)+R01(2,14)+R01(2,18)
         F(1,2,3,3)= R04( 9)-R03(5,1)-R03(5,4)+R03(9,7)*QX
     *              +R02(4,4)+R02(4,5)-R02(6,16)*QX-R02(6,20)*QX
     *              +R01(2,15)*QX+R01(2,19)*QX
         F(2,2,3,3)= R04(13)-R03(8,1)-R03(8,4)+R02(2,4)+R02(2,5)
     *              +R02(3,26)-R01(3,35)-R01(3,40)+R00(20)+R00(25)
         F(3,2,3,3)= R04(14)-R03(9,1)-R03(9,4)+R03(9,7)*QZ
     *              +R02(6,4)+R02(6,5)-R02(6,16)*QZ-R02(6,20)*QZ
     *              +R01(2,15)*QZ+R01(2,19)*QZ
         F(0,3,3,3)= R03(10,6)-R02(3,15)-R02(3,19)+R02(3,24)*QZ
     *              +R01(3,14)+R01(3,18)-R01(3,33)*QZ-R01(3,38)*QZ
     *              +R00(18)*QZ+R00(23)*QZ
         F(1,3,3,3)= R04(10)-R03(6,1)-R03(6,4)+R03(10,7)*QX+R03(6,8)*QZ
     *              +R02(5,4)+R02(5,5)-R02(3,16)*QX-R02(5,17)*QZ
     *              -R02(3,20)*QX-R02(5,21)*QZ+R02(3,25)*QXZ
     *              +R01(3,15)*QX+R01(1,16)*QZ+R01(3,19)*QX+R01(1,20)*QZ
     *              -R01(3,34)*QXZ-R01(3,39)*QXZ+R00(19)*QXZ+R00(24)*QXZ
         F(2,3,3,3)= R04(14)-R03(9,1)-R03(9,4)+R03(9,8)*QZ
     *              +R02(6,4)+R02(6,5)-R02(6,17)*QZ-R02(6,21)*QZ
     *              +R01(2,16)*QZ+R01(2,20)*QZ
         F(3,3,3,3)=R04(15)-R03(10,1)-R03(10,4)+(R03(10,7)+R03(10,8))*QZ
     *              +R02(3,4)+R02(3,5)-R02(3,16)*QZ-R02(3,17)*QZ
     *              -R02(3,20)*QZ-R02(3,21)*QZ+R02(3,25)*QZZ+R02(3,26)
     *              +R01(3,15)*QZ+R01(3,16)*QZ+R01(3,19)*QZ+R01(3,20)*QZ
     *              -R01(3,34)*QZZ-R01(3,35)-R01(3,39)*QZZ-R01(3,40)
     *              +R00(19)*QZZ+R00(20)+R00(24)*QZZ+R00(25)
      ENDIF
      RETURN
      END
C
C*MODULE INT2B   *DECK R30S1S
      SUBROUTINE R30S1S(JTYPE,F,P)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C ROTATE UP TO 256 INTEGRALS TO SPACE FIXED AXES
C INCOMING AND OUTGOING INTEGRALS IN F
C P(1,1),...ARE DIRECTION COSINES OF SPACE FIXED AXES WRT AXES AT P
C
      DIMENSION F(0:3,0:3,0:3,0:3),P(3,3)
C
      DIMENSION T(3)
C
      IF(JTYPE.EQ.6) THEN
         DO 610 K=0,3
            DO 610 J=0,3
               DO 610 I=0,3
                  T(1)= F(I,J,K,1)
                  T(2)= F(I,J,K,2)
                  T(3)= F(I,J,K,3)
                  F(I,J,K,1)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
                  F(I,J,K,2)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
                  F(I,J,K,3)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  610    CONTINUE
         DO 620 L=0,3
            DO 620 J=0,3
               DO 620 I=0,3
                  T(1)= F(I,J,1,L)
                  T(2)= F(I,J,2,L)
                  T(3)= F(I,J,3,L)
                  F(I,J,1,L)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
                  F(I,J,2,L)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
                  F(I,J,3,L)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  620    CONTINUE
         DO 630 L=0,3
            DO 630 K=0,3
               DO 630 I=0,3
                  T(1)= F(I,1,K,L)
                  T(2)= F(I,2,K,L)
                  T(3)= F(I,3,K,L)
                  F(I,1,K,L)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
                  F(I,2,K,L)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
                  F(I,3,K,L)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  630    CONTINUE
         DO 640 L=0,3
            DO 640 K=0,3
               DO 640 J=0,3
                  T(1)= F(1,J,K,L)
                  T(2)= F(2,J,K,L)
                  T(3)= F(3,J,K,L)
                  F(1,J,K,L)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
                  F(2,J,K,L)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
                  F(3,J,K,L)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  640    CONTINUE
      ELSEIF(JTYPE.EQ.5) THEN
         DO 510 J=0,3
            DO 510 I=0,3
               T(1)= F(I,J,1,0)
               T(2)= F(I,J,2,0)
               T(3)= F(I,J,3,0)
               F(I,J,1,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
               F(I,J,2,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
               F(I,J,3,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  510    CONTINUE
         DO 520 K=0,3
            DO 520 I=0,3
               T(1)= F(I,1,K,0)
               T(2)= F(I,2,K,0)
               T(3)= F(I,3,K,0)
               F(I,1,K,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
               F(I,2,K,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
               F(I,3,K,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  520    CONTINUE
         DO 530 K=0,3
            DO 530 J=0,3
               T(1)= F(1,J,K,0)
               T(2)= F(2,J,K,0)
               T(3)= F(3,J,K,0)
               F(1,J,K,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
               F(2,J,K,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
               F(3,J,K,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  530    CONTINUE
      ELSEIF(JTYPE.EQ.4) THEN
         DO 410 I=0,3
            T(1)= F(I,0,1,0)
            T(2)= F(I,0,2,0)
            T(3)= F(I,0,3,0)
            F(I,0,1,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
            F(I,0,2,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
            F(I,0,3,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  410    CONTINUE
         DO 420 K=0,3
            T(1)= F(1,0,K,0)
            T(2)= F(2,0,K,0)
            T(3)= F(3,0,K,0)
            F(1,0,K,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
            F(2,0,K,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
            F(3,0,K,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  420    CONTINUE
      ELSEIF(JTYPE.EQ.3) THEN
         DO 310 I=0,3
            T(1)= F(I,1,0,0)
            T(2)= F(I,2,0,0)
            T(3)= F(I,3,0,0)
            F(I,1,0,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
            F(I,2,0,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
            F(I,3,0,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  310    CONTINUE
         DO 320 J=0,3
            T(1)= F(1,J,0,0)
            T(2)= F(2,J,0,0)
            T(3)= F(3,J,0,0)
            F(1,J,0,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
            F(2,J,0,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
            F(3,J,0,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
  320    CONTINUE
      ELSEIF(JTYPE.EQ.2) THEN
         T(1)= F(1,0,0,0)
         T(2)= F(2,0,0,0)
         T(3)= F(3,0,0,0)
         F(1,0,0,0)= T(1)*P(1,1)+T(2)*P(2,1)+T(3)*P(3,1)
         F(2,0,0,0)= T(1)*P(1,2)+T(2)*P(2,2)+T(3)*P(3,2)
         F(3,0,0,0)= T(1)*P(1,3)+T(2)*P(2,3)+T(3)*P(3,3)
      ENDIF
      RETURN
      END
C
C*MODULE INT2B   *DECK SP0S1S
      SUBROUTINE SP0S1S(JTYPE)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (ZER=0.0D+00, PT5=0.5D+00)
C
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
C
C     ==================================================================
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(2,9),FQ2(27),FQ3(20),FQ4(5)
      COMMON /JMSGYH/ SQ(4)
      COMMON /KI2   / ACY,ACY2,AQX,AQX2,AQXY,Y03,Y04
      COMMON /KI3   / R00(25),R01(3,40),R02(156),R03(80),R04(15)
C
      DIMENSION WORK(4)
C
C FROM HERE, MCMURCHIE-DAVIDSON ALGORITHM IS USED.
C
C GENERATE (R] INTEGRALS
C
      IF(JTYPE.EQ.6) THEN
         DO I=1,5
            FQ0(I)= ZER
         ENDDO
         DO  J=1,9
               FQ1(1,J)= ZER
               FQ1(2,J)= ZER
         ENDDO
         DO  J=1,27
               FQ2(J)= ZER
         ENDDO
         DO  I=1,20
               FQ3(I)= ZER
         ENDDO
         DO  I=1,5
            FQ4(I)= ZER
         ENDDO
         CALL INTJ6
         CALL INTK6
C
      ELSEIF(JTYPE.EQ.5) THEN
         FQ0(1)= ZER
         FQ0(2)= ZER
         DO 150 J=1,3
               FQ1(1,J)= ZER
               FQ1(2,J)= ZER
  150    CONTINUE
         DO 160 J=1,9
               FQ2(J)= ZER
  160    CONTINUE
            FQ3(1)= ZER
            FQ3(2)= ZER
            FQ3(3)= ZER
            FQ3(4)= ZER
         CALL INTJ5
         CALL INTK5
C
      ELSEIF(JTYPE.EQ.4) THEN
         FQ0(1)= ZER
         FQ0(2)= ZER
         DO 120 J=1,3
               FQ1(1,J)= ZER
               FQ1(2,J)= ZER
  120    CONTINUE
         FQ2(1)= ZER
         FQ2(2)= ZER
         FQ2(3)= ZER
         CALL INTJ4
C
      XMD2= X43*PT5
         WORK(1)= XMD2*SQ(2)
         WORK(2)= SQ(1)
         WORK(3)= WORK(1)
         WORK(4)=-Y03*SQ(2)
         R00(1)= R00(1)+FQ0(1)*WORK(2)
         R00(2)= R00(2)+FQ0(1)*WORK(4)
         R00(3)= R00(3)+FQ0(2)*WORK(2)
         R00(4)= R00(4)+FQ0(2)*WORK(4)
C
         DO I=1,3
            R01(1,I)= R01(1,I)-FQ1(1,I)*AQX*WORK(I)
            R01(2,I)= R01(2,I)-FQ1(1,I)*ACY*WORK(I)
            R01(3,I)= R01(3,I)+FQ1(2,I)*WORK(I)
         ENDDO
         R01(1,4)= R01(1,4)-FQ1(1,2)*AQX*WORK(4)
         R01(2,4)= R01(2,4)-FQ1(1,2)*ACY*WORK(4)
         R01(3,4)= R01(3,4)+FQ1(2,2)*WORK(4)
C
         R02(1)= R02(1)+(FQ2(1)*AQX2-FQ1(1,2))*WORK(1)
         R02(2)= R02(2)+(FQ2(1)*ACY2-FQ1(1,2))*WORK(1)
         R02(3)= R02(3)+(FQ2(3)     -FQ1(1,2))*WORK(1)
         R02(4)= R02(4)+FQ2(1)*AQXY*WORK(1)
         R02(5)= R02(5)-FQ2(2)*AQX *WORK(1)
         R02(6)= R02(6)-FQ2(2)*ACY *WORK(1)
      ELSEIF(JTYPE.EQ.3) THEN
         FQ0(1)= ZER
         FQ1(1,1)= ZER
         FQ1(2,1)= ZER
         FQ2(1)= ZER
         FQ2(2)= ZER
         FQ2(3)= ZER
         CALL INTJ3
C
      XMD2= X43*PT5
         R00(1)= R00(1)+FQ0(1)*SQ(1)
         R00(2)= R00(2)-FQ0(1)*Y03*SQ(2)
         R00(3)= R00(3)+FQ0(1)*Y04*SQ(3)
         R00(4)= R00(4)+FQ0(1)*XMD2*SQ(4)
         R00(5)= R00(5)-FQ0(1)*Y03*Y04*SQ(4)
C
         WORK(1)= XMD2*SQ(2)
         WORK(2)= XMD2*SQ(3)
         WORK(3)=-Y03*XMD2*SQ(4)
         WORK(4)= Y04*XMD2*SQ(4)
C        FQ11=-FQ1(1,1)*AQX
C        FQ12=-FQ1(1,1)*ACY
C        FQ13= FQ1(2,1)
         DO I=1,4
C           R01(1,I)= R01(1,I)+FQ11*WORK(I)
C           R01(2,I)= R01(2,I)+FQ12*WORK(I)
C           R01(3,I)= R01(3,I)+FQ13*WORK(I)
            R01(1,I)= R01(1,I)-FQ1(1,1)*AQX*WORK(I)
            R01(2,I)= R01(2,I)-FQ1(1,1)*ACY*WORK(I)
            R01(3,I)= R01(3,I)+FQ1(2,1)*WORK(I)
         ENDDO
C
         XMD4= XMD2*XMD2*SQ(4)
         R02(1)= R02(1)+(FQ2(1)*AQX2-FQ1(1,1))*XMD4
         R02(2)= R02(2)+(FQ2(1)*ACY2-FQ1(1,1))*XMD4
         R02(3)= R02(3)+(FQ2(3)-FQ1(1,1))*XMD4
         R02(4)= R02(4)+FQ2(1)*AQXY*XMD4
         R02(5)= R02(5)-FQ2(2)*AQX*XMD4
         R02(6)= R02(6)-FQ2(2)*ACY*XMD4
      ELSEIF(JTYPE.EQ.2) THEN
         FQ0(1)= ZER
         FQ1(1,1)= ZER
         FQ1(2,1)= ZER
         CALL INTJ2
C
      XMD2= X43*PT5*SQ(2)
         R00(1)= R00(1)+FQ0(1)*SQ(1)
         R00(2)= R00(2)-FQ0(1)*Y03*SQ(2)
         R01(1,1)= R01(1,1)-FQ1(1,1)*AQX*XMD2
         R01(2,1)= R01(2,1)-FQ1(1,1)*ACY*XMD2
         R01(3,1)= R01(3,1)+FQ1(2,1)*XMD2
      ELSEIF(JTYPE.EQ.1) THEN
         FQ0(1)= ZER
         CALL INTJ1
C
         R00(1)= R00(1)+FQ0(1)*SQ(1)
      ENDIF
      RETURN
      END
C
C*MODULE INT2B   *DECK INTJ1
      SUBROUTINE INTJ1
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=1 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (ONE=1.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(18),FQ2(27),FQ3(20),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D00P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         XVA=(PQS+QPS)*X12*X34*X41
         IF(LRINT) THEN
            RHO0=X12*X34*X41
            RHO =RHO0*EMU2/(RHO0+EMU2)
            XVA =(PQS+QPS)*RHO
         END IF
C
C            FM(T) EVALUATION
C
         IF (XVA.LE.TMAX) THEN
            TVAL=XVA*RFINC(0)
            IPT=NINT(TVAL)
            FX =      FGRID(4,IPT,0) *TVAL
            FX = (FX +FGRID(3,IPT,0))*TVAL
            FX = (FX +FGRID(2,IPT,0))*TVAL
            FX = (FX +FGRID(1,IPT,0))*TVAL
            FX =  FX +FGRID(0,IPT,0)
            FQ(0)=FX
            FQ(0) = FQ(0)*FQ00*SQRT(X41)
         ELSE
            XIN= PI4*X41/XVA
            FQ(0)= SQRT(XIN)*FQ00
         ENDIF
         IF(LRINT)FQ(0)=FQ(0)*SQRT(RHO/RHO0)
         FQ0(1)= FQ0(1)+FQ(0)
  300 CONTINUE
      RETURN
      END
C
C*MODULE INT2B   *DECK INTJ2
      SUBROUTINE INTJ2
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=2 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(18),FQ2(27),FQ3(20),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D00P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         RHO= X12*X34*X41
         IF(LRINT) THEN
            RHO0 = RHO
            RHO  = RHO0*EMU2/(RHO0+EMU2)
         ENDIF
         XVA=(PQS+QPS)*RHO
C
C            FM(T) EVALUATION
C
         IF (XVA.LE.TMAX) THEN
            TVAL=XVA*RFINC(0)
            IPT=NINT(TVAL)
            FX =      FGRID(4,IPT,0) *TVAL
            FX = (FX +FGRID(3,IPT,0))*TVAL
            FX = (FX +FGRID(2,IPT,0))*TVAL
            FX = (FX +FGRID(1,IPT,0))*TVAL
            FX =  FX +FGRID(0,IPT,0)
            FQ(0)=FX
            TVAL=XVA*RFINC(1)
            IPT=NINT(TVAL)
            FX =      FGRID(4,IPT,1) *TVAL
            FX = (FX +FGRID(3,IPT,1))*TVAL
            FX = (FX +FGRID(2,IPT,1))*TVAL
            FX = (FX +FGRID(1,IPT,1))*TVAL
            FX =  FX +FGRID(0,IPT,1)
            FQ(1)=FX
            FQS = FQ00*SQRT(X41)
            FQ(0) = FQ(0)*FQS
            FQ(1) = FQ(1)*FQS*RHO*TWO
         ELSE
            XIN= ONE/XVA
            FQ(0)= SQRT(PI4*XIN*X41)*FQ00
            FQ(1)= FQ(0)*XIN*RHO
         ENDIF
         IF(LRINT) THEN
            FQ(0)=FQ(0)*SQRT(RHO/RHO0)
            FQ(1)=FQ(1)*SQRT(RHO/RHO0)
         ENDIF
C
         FQ0(1)= FQ0(1)+FQ(0)
         FQ1(1)= FQ1(1)+FQ(1)
         FQ1(2)= FQ1(2)+FQ(1)*PQR
  300 CONTINUE
      RETURN
      END
C
C*MODULE INT2B   *DECK INTJ3
      SUBROUTINE INTJ3
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=3 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(18),FQ2(27),FQ3(20),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D00P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         RHO= X12*X34*X41
         IF(LRINT) THEN
            RHO0 = RHO
            RHO  = RHO0*EMU2/(RHO0+EMU2)
         ENDIF
         XVA=(PQS+QPS)*RHO
         RHO= RHO*TWO
C
C            FM(T) EVALUATION
C
         IF (XVA.LE.TMAX) THEN
           TVAL=XVA*RFINC(2)
           IPT=NINT(TVAL)
           FX =      FGRID(4,IPT,2) *TVAL
           FX = (FX +FGRID(3,IPT,2))*TVAL
           FX = (FX +FGRID(2,IPT,2))*TVAL
           FX = (FX +FGRID(1,IPT,2))*TVAL
           FX =  FX +FGRID(0,IPT,2)
           FQ(2)=FX
           TVAL=XVA*RXINC
           IPT=NINT(TVAL)
           ET =       XGRID(4,IPT) *TVAL
           ET = (ET + XGRID(3,IPT))*TVAL
           ET = (ET + XGRID(2,IPT))*TVAL
           ET = (ET + XGRID(1,IPT))*TVAL
           ET =  ET + XGRID(0,IPT)
           T2 = XVA*TWO
           FQ( 1)=RMR( 2)*(T2*FQ( 2)+ET)
           FQ( 0)=RMR( 1)*(T2*FQ( 1)+ET)
           FQS = FQ00*SQRT(X41)
           IF(LRINT) FQS=FQS*SQRT(RHO/RHO0/TWO)
           DO MXX = 0, 2
             FQ(MXX) = FQ(MXX)*FQS
             FQS = FQS*RHO
           END DO
         ELSE
            XIN= ONE/XVA
            FQ(0)= SQRT(PI4*XIN*X41)*FQ00
            IF(LRINT)FQ(0)= FQ(0)*SQRT(RHO/RHO0/TWO)
            TMP= PT5*XIN
            DO 220 M=1,2
               FQ(M)= FQ(M-1)*TMP*RHO
  220       TMP= TMP+XIN
         ENDIF
C
         FQ0(1)= FQ0(1)+FQ(0)
         FQ1(1)= FQ1(1)+FQ(1)
         FQ1(2)= FQ1(2)+FQ(1)*PQR
         FQ2(1)= FQ2(1)+FQ(2)
         FQ2(2)= FQ2(2)+FQ(2)*PQR
         FQ2(3)= FQ2(3)+FQ(2)*PQS
  300 CONTINUE
      RETURN
      END
C*MODULE INT2B   *DECK INTJ4
      SUBROUTINE INTJ4
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=4 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(18),FQ2(27),FQ3(20),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D01P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         RHO= X12*X34*X41
         IF(LRINT) THEN
           RHO0 = RHO
           RHO  = RHO0*EMU2/(RHO0+EMU2)
         ENDIF
         XVA=(PQS+QPS)*RHO
         RHO= RHO*TWO
C
C            FM(T) EVALUATION
C
         IF (XVA.LE.TMAX) THEN
           TVAL=XVA*RFINC(2)
           IPT=NINT(TVAL)
           FX =      FGRID(4,IPT,2) *TVAL
           FX = (FX +FGRID(3,IPT,2))*TVAL
           FX = (FX +FGRID(2,IPT,2))*TVAL
           FX = (FX +FGRID(1,IPT,2))*TVAL
           FX =  FX +FGRID(0,IPT,2)
           FQ(2)=FX
           TVAL=XVA*RXINC
           IPT=NINT(TVAL)
           ET =       XGRID(4,IPT) *TVAL
           ET = (ET + XGRID(3,IPT))*TVAL
           ET = (ET + XGRID(2,IPT))*TVAL
           ET = (ET + XGRID(1,IPT))*TVAL
           ET =  ET + XGRID(0,IPT)
           T2 = XVA*TWO
           FQ( 1)=RMR( 2)*(T2*FQ( 2)+ET)
           FQ( 0)=RMR( 1)*(T2*FQ( 1)+ET)
           FQS = FQ00*SQRT(X41)
            IF(LRINT) FQS=FQS*SQRT(RHO/RHO0/TWO)
           DO MXX = 0, 2
             FQ(MXX) = FQ(MXX)*FQS
             FQS = FQS*RHO
           END DO
         ELSE
            XIN= ONE/XVA
            FQ(0)= SQRT(PI4*XIN*X41)*FQ00
            IF(LRINT)FQ(0)= FQ(0)*SQRT(RHO/RHO0/TWO)
            TMP= PT5*XIN
            DO 220 M=1,2
               FQ(M)= FQ(M-1)*TMP*RHO
  220       TMP= TMP+XIN
         ENDIF
C
         XMD1= TX21(I)
         Y01= TY01(I)
         DP00= D00P(I)
         TMP1= XMD1*PQR
         TMP2= XMD1*PQS
         FQ0(1)= FQ0(1)+FQ(0)*DP00
         FQ0(2)= FQ0(2)+FQ(0)*Y01
         FQ1(1)= FQ1(1)+FQ(1)*DP00
         FQ1(2)= FQ1(2)+FQ(1)*DP00*PQR
         FQ1(3)= FQ1(3)+FQ(1)*XMD1
         FQ1(4)= FQ1(4)+FQ(1)*TMP1
         FQ1(5)= FQ1(5)+FQ(1)*Y01
         FQ1(6)= FQ1(6)+FQ(1)*Y01*PQR
         FQ2(1)= FQ2(1)+FQ(2)*XMD1
         FQ2(2)= FQ2(2)+FQ(2)*TMP1
         FQ2(3)= FQ2(3)+FQ(2)*TMP2
  300 CONTINUE
      RETURN
      END
C
C*MODULE INT2B   *DECK INTJ5
      SUBROUTINE INTJ5
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=5 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER ( PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(2,9),FQ2(3,9),FQ3(20),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DIMENSION WORK(3,3)
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D01P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         RHO= X12*X34*X41
         IF(LRINT) THEN
           RHO0 = RHO
           RHO  = RHO0*EMU2/(RHO0+EMU2)
         ENDIF
         XVA=(PQS+QPS)*RHO
         RHO= RHO*TWO
C
C            FM(T) EVALUATION...DOWNWARD RECURSION FOR JTYPE.GE.5
C
         IF (XVA.LE.TMAX) THEN
           TVAL=XVA*RFINC(3)
           IPT=NINT(TVAL)
           FX =      FGRID(4,IPT,3) *TVAL
           FX = (FX +FGRID(3,IPT,3))*TVAL
           FX = (FX +FGRID(2,IPT,3))*TVAL
           FX = (FX +FGRID(1,IPT,3))*TVAL
           FX =  FX +FGRID(0,IPT,3)
           FQ(3)=FX
           TVAL=XVA*RXINC
           IPT=NINT(TVAL)
           ET =      XGRID(4,IPT) *TVAL
           ET = (ET +XGRID(3,IPT))*TVAL
           ET = (ET +XGRID(2,IPT))*TVAL
           ET = (ET +XGRID(1,IPT))*TVAL
           ET =  ET +XGRID(0,IPT)
           T2 = XVA*TWO
           DO MXX = 3, 1, -1
             FQ(MXX-1) = RMR(MXX)*(T2*FQ(MXX)+ET)
           END DO
           FQS = FQ00*SQRT(X41)
            IF(LRINT) FQS=FQS*SQRT(RHO/RHO0/TWO)
           DO MXX = 0, 3
             FQ(MXX) = FQ(MXX)*FQS
             FQS = FQS*RHO
           END DO
         ELSE
            XIN= ONE/XVA
            FQ(0)= SQRT(PI4*XIN*X41)*FQ00
            IF(LRINT)FQ(0)= FQ(0)*SQRT(RHO/RHO0/TWO)
            TMP= PT5*XIN
            DO 220 M=1,3
               FQ(M)= FQ(M-1)*TMP*RHO
  220       TMP= TMP+XIN
         ENDIF
C
         XMD1= TX21(I)
         Y01= TY01(I)
         DP00= D00P(I)
         WORK(1,1)= DP00
         WORK(1,2)= Y01
         WORK(1,3)= XMD1
         DO J=1,3
          WORK(2,J)= WORK(1,J)*PQR
          WORK(3,J)= WORK(1,J)*PQS
         ENDDO
C
         FQ0(1)= FQ0(1)+FQ(0)*WORK(1,1)
         FQ0(2)= FQ0(2)+FQ(0)*WORK(1,2)
         DO J=1,3
          FQ1(1,J)= FQ1(1,J)+FQ(1)*WORK(1,J)
          FQ1(2,J)= FQ1(2,J)+FQ(1)*WORK(2,J)
         ENDDO
         DO 230 J=1,3
            FQ2(1,J)= FQ2(1,J)+FQ(2)*WORK(1,J)
            FQ2(2,J)= FQ2(2,J)+FQ(2)*WORK(2,J)
            FQ2(3,J)= FQ2(3,J)+FQ(2)*WORK(3,J)
  230    CONTINUE
            FQ3(1)= FQ3(1)+FQ(3)*WORK(1,3)
            FQ3(2)= FQ3(2)+FQ(3)*WORK(2,3)
            FQ3(3)= FQ3(3)+FQ(3)*WORK(3,3)
            FQ3(4)= FQ3(4)+FQ(3)*WORK(3,3)*PQR
  300 CONTINUE
      RETURN
      END
C
C*MODULE INT2B   *DECK INTK5
      SUBROUTINE INTK5
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=5 INTEGRALS
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (PT5=0.5D+00, THREE=3.0D+00)
C
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(2,9),FQ2(3,9),FQ3(4,5),FQ4(5)
      COMMON /JMSGYH/ SQ(4)
      COMMON /KI2   / ACY,ACY2,AQX,AQX2,AQXY,Y03,Y04
      COMMON /KI3   / R00(5,5),R01(3,40),R02(6,26),R03(10,8),R04(15)
C
      DIMENSION WORK(9),WORK2(6,3)
         XMD2= X43*PT5
         WORK(1)= SQ(1)
         WORK(2)=-Y03*SQ(2)
         WORK(3)= Y04*SQ(3)
         WORK(4)=-Y03*Y04*SQ(4)
         WORK(5)= XMD2*SQ(4)
         WORK(6)= XMD2*SQ(2)
         WORK(7)= XMD2*SQ(3)
         WORK(8)=-XMD2*Y03*SQ(4)
         WORK(9)= XMD2*Y04*SQ(4)
         XMD4= XMD2*XMD2*SQ(4)
         DO I=1,2
          DO J=1,5
            R00( J,I)= R00( J,I)+FQ0(I)*WORK(J)
          ENDDO
         ENDDO
C
         DO I=1,3
            WORK2(1,I)=-FQ1(1,I)*AQX
            WORK2(2,I)=-FQ1(1,I)*ACY
            WORK2(3,I)= FQ1(2,I)
         ENDDO
         DO I=1,4
            R01(1,I)= R01(1,I)+WORK2(1,1)*WORK(I+5)
            R01(2,I)= R01(2,I)+WORK2(2,1)*WORK(I+5)
            R01(3,I)= R01(3,I)+WORK2(3,1)*WORK(I+5)
         ENDDO
         DO I=5,8
            R01(1,I)= R01(1,I)+WORK2(1,2)*WORK(I+1)
            R01(2,I)= R01(2,I)+WORK2(2,2)*WORK(I+1)
            R01(3,I)= R01(3,I)+WORK2(3,2)*WORK(I+1)
         ENDDO
         DO I=9,13
            R01(1,I)= R01(1,I)+WORK2(1,3)*WORK(I-8)
            R01(2,I)= R01(2,I)+WORK2(2,3)*WORK(I-8)
            R01(3,I)= R01(3,I)+WORK2(3,3)*WORK(I-8)
         ENDDO
C
         XMD4XY=XMD4*AQXY
         XMD4X=-XMD4*AQX
         XMD4Y=-XMD4*ACY
         DO I=1,2
            R02(1,I)= R02(1,I)+(FQ2(1,I)*AQX2-FQ1(1,I))*XMD4
            R02(2,I)= R02(2,I)+(FQ2(1,I)*ACY2-FQ1(1,I))*XMD4
            R02(3,I)= R02(3,I)+(FQ2(3,I)     -FQ1(1,I))*XMD4
            R02(4,I)= R02(4,I)+FQ2(1,I)*XMD4XY
            R02(5,I)= R02(5,I)+FQ2(2,I)*XMD4X
            R02(6,I)= R02(6,I)+FQ2(2,I)*XMD4Y
         ENDDO
         WORK2(1,1)= FQ2(1,3)*AQX2-FQ1(1,3)
         WORK2(2,1)= FQ2(1,3)*ACY2-FQ1(1,3)
         WORK2(3,1)= FQ2(3,3)     -FQ1(1,3)
         WORK2(4,1)= FQ2(1,3)*AQXY
         WORK2(5,1)=-FQ2(2,3)*AQX
         WORK2(6,1)=-FQ2(2,3)*ACY
         DO I=3,6
          DO J=1,6
            R02(J,I)= R02(J,I)+WORK2(J,1)*WORK(I+3)
          ENDDO
         ENDDO
C
         R03( 1,1)= R03( 1,1)+(FQ3(1,1)*AQX2-FQ2(1,3)*THREE)*XMD4X
         R03( 2,1)= R03( 2,1)+(FQ3(1,1)*AQX2-FQ2(1,3)      )*XMD4Y
         R03( 3,1)= R03( 3,1)+(FQ3(2,1)*AQX2-FQ2(2,3)      )*XMD4
         R03( 4,1)= R03( 4,1)+(FQ3(1,1)*ACY2-FQ2(1,3)      )*XMD4X
         R03( 5,1)= R03( 5,1)+(FQ3(2,1)*AQXY               )*XMD4
         R03( 6,1)= R03( 6,1)+(FQ3(3,1)     -FQ2(1,3)      )*XMD4X
         R03( 7,1)= R03( 7,1)+(FQ3(1,1)*ACY2-FQ2(1,3)*THREE)*XMD4Y
         R03( 8,1)= R03( 8,1)+(FQ3(2,1)*ACY2-FQ2(2,3)      )*XMD4
         R03( 9,1)= R03( 9,1)+(FQ3(3,1)     -FQ2(1,3)      )*XMD4Y
         R03(10,1)= R03(10,1)+(FQ3(4,1)     -FQ2(2,3)*THREE)*XMD4
      RETURN
      END
C
C*MODULE INT2B   *DECK INTJ6
      SUBROUTINE INTJ6
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=6 INTEGRALS
C
      LOGICAL   LCFLAG,LRINT
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER ( PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI4=0.785398163397448D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(2,9),FQ2(3,9),FQ3(4,5),FQ4(5)
      PARAMETER (NTX=4, NPF=450, NGRD=7, NPX=1000, MXQT=16)
      COMMON /FMTTBL/ FGRID(0:NTX,0:NPF,0:NGRD),XGRID(0:NTX,0:NPX)
     *,               TMAX,RFINC(0:NGRD),RXINC
     *,               RMR(MXQT),TLGM(0:MXQT),NORD
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /MAXC  / CMAX(MXGTOT),CMAXA(MXGSH),CMAXB(MXGSH),
     2                CMAXC(MXGSH),CMAXD(MXGSH),ISMLP(MXG2),ISMLQ
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
C
      DIMENSION WORK(4,9)
C
      DO 300 I=1,NGANGB
         ISML= ISMLQ+ISMLP(I)
         IF(ISML.GE.2) GO TO 300
         X12= TX12(I)
         Y02= TY02(I)
         FQ00= D11P(I)
         X41= ONE/(X12+X34)
         PQR= Y02-AQZ
         PQS= PQR*PQR
         RHO= X12*X34*X41
         IF(LRINT) THEN
            RHO0 = RHO
            RHO  = RHO0*EMU2/(RHO0+EMU2)
         ENDIF
         XVA=(PQS+QPS)*RHO
         RHO= RHO*TWO
C
C            FM(T) EVALUATION...DOWNWARD RECURSION FOR JTYPE.GE.5
C
         IF (XVA.LE.TMAX) THEN
           TVAL=XVA*RFINC(4)
           IPT=NINT(TVAL)
           FX =      FGRID(4,IPT,4) *TVAL
           FX = (FX +FGRID(3,IPT,4))*TVAL
           FX = (FX +FGRID(2,IPT,4))*TVAL
           FX = (FX +FGRID(1,IPT,4))*TVAL
           FX =  FX +FGRID(0,IPT,4)
           FQ(4)=FX
           TVAL=XVA*RXINC
           IPT=NINT(TVAL)
           ET =      XGRID(4,IPT) *TVAL
           ET = (ET +XGRID(3,IPT))*TVAL
           ET = (ET +XGRID(2,IPT))*TVAL
           ET = (ET +XGRID(1,IPT))*TVAL
           ET =  ET +XGRID(0,IPT)
           T2 = XVA*TWO
           DO MXX = 4, 1, -1
             FQ(MXX-1) = RMR(MXX)*(T2*FQ(MXX)+ET)
           END DO
           FQS = FQ00*SQRT(X41)
            IF(LRINT) FQS=FQS*SQRT(RHO/RHO0/TWO)
           DO MXX = 0, 4
             FQ(MXX) = FQ(MXX)*FQS
             FQS = FQS*RHO
           END DO
         ELSE
            XIN= ONE/XVA
            FQ(0)= SQRT(PI4*XIN*X41)*FQ00
            IF(LRINT)FQ(0)= FQ(0)*SQRT(RHO/RHO0/TWO)
            TMP= PT5*XIN
            DO 220 M=1,4
               FQ(M)= FQ(M-1)*TMP*RHO
  220       TMP= TMP+XIN
         ENDIF
C
         XMD1= TX21(I)
         Y01= TY01(I)
         DP00= D00P(I)
         DP01= D01P(I)
         DP10= D10P(I)
         PQT = PQR*PQS
         WORK(1,1)= DP00
         WORK(1,2)= Y01*DP01
         WORK(1,3)= Y02*DP10
         WORK(1,4)= Y01*Y02
         WORK(1,5)= XMD1
         WORK(1,6)= XMD1*DP01
         WORK(1,7)= XMD1*DP10
         WORK(1,8)= XMD1*Y01
         WORK(1,9)= XMD1*XMD1
         DO 230 J=1,4
            WORK(2,J)= WORK(1,J)*PQR
            WORK(3,J)= WORK(1,J)*PQS
  230    CONTINUE
         DO 240 J=5,9
            WORK(2,J)= WORK(1,J)*PQR
            WORK(3,J)= WORK(1,J)*PQS
            WORK(4,J)= WORK(1,J)*PQT
  240    CONTINUE
C
         DO 250 J=1,5
            FQ0(J)= FQ0(J)+FQ(0)*WORK(1,J)
  250    CONTINUE
C
         DO 260 J=1,8
            FQ1(1,J)= FQ1(1,J)+FQ(1)*WORK(1,J)
            FQ1(2,J)= FQ1(2,J)+FQ(1)*WORK(2,J)
  260    CONTINUE
         FQ1(1,9)= FQ1(1,9)+FQ(1)*WORK(1,9)
C
         DO 270 J=1,9
C           DO 270 K=1,3
C              FQ2(K,J)= FQ2(K,J)+FQ(2)*WORK(K,J)
               FQ2(1,J)= FQ2(1,J)+FQ(2)*WORK(1,J)
               FQ2(2,J)= FQ2(2,J)+FQ(2)*WORK(2,J)
               FQ2(3,J)= FQ2(3,J)+FQ(2)*WORK(3,J)
  270    CONTINUE
C
         DO 280 J=1,5
C           DO 280 K=1,4
C              FQ3(K,J)= FQ3(K,J)+FQ(3)*WORK(K,J+4)
               FQ3(1,J)= FQ3(1,J)+FQ(3)*WORK(1,J+4)
               FQ3(2,J)= FQ3(2,J)+FQ(3)*WORK(2,J+4)
               FQ3(3,J)= FQ3(3,J)+FQ(3)*WORK(3,J+4)
               FQ3(4,J)= FQ3(4,J)+FQ(3)*WORK(4,J+4)
  280    CONTINUE
C
         FQ4(1)= FQ4(1)+FQ(4)*WORK(1,9)
         FQ4(2)= FQ4(2)+FQ(4)*WORK(2,9)
         FQ4(3)= FQ4(3)+FQ(4)*WORK(3,9)
         FQ4(4)= FQ4(4)+FQ(4)*WORK(4,9)
         FQ4(5)= FQ4(5)+FQ(4)*WORK(4,9)*PQR
  300 CONTINUE
      RETURN
      END
C*MODULE INT2B   *DECK INTK6
      SUBROUTINE INTK6
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C GENERATE JTYPE=6 INTEGRALS
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      PARAMETER (PT5=0.5D+00, THREE=3.0D+00, SIX=6.0D+00)
C
      COMMON /FQ04  / FQ(0:4),FQ0(5),FQ1(2,9),FQ2(3,9),FQ3(4,5),FQ4(5)
      COMMON /GEOMPQ/ R12,RAB,X34,X43,AQZ,QPR,QPS,
     2                TX12(MXG2),TX21(MXG2),TY01(MXG2),TY02(MXG2),
     3                D00P(MXG2),D01P(MXG2),D10P(MXG2),D11P(MXG2),
     4                NGANGB
      COMMON /JMSGYH/ SQ(4)
      COMMON /KI2   / ACY,ACY2,AQX,AQX2,AQXY,Y03,Y04
      COMMON /KI3   / R00(5,5),R01(3,40),R02(6,26),R03(10,8),R04(15)
C
      DIMENSION  WORK(9),FQWK(6,9),WORK2(10)
C
      XMD2= X43*PT5
      XMD4= XMD2*XMD2*SQ(4)
         WORK( 1)= SQ(1)
         WORK( 2)=-Y03*SQ(2)
         WORK( 3)= Y04*SQ(3)
         WORK( 4)=-Y03*Y04*SQ(4)
         WORK( 5)= XMD2*SQ(4)
         WORK( 6)= XMD2*SQ(2)
         WORK( 7)= XMD2*SQ(3)
         WORK( 8)=-XMD2*Y03*SQ(4)
         WORK( 9)= XMD2*Y04*SQ(4)
C
C
      DO 10 J=1,5
         DO 10 I=1,5
            R00(I,J)= R00(I,J)+FQ0(J)*WORK(I)
   10 CONTINUE
C
      DO I=1,8
       FQWK(1,I)=-FQ1(1,I)*AQX
       FQWK(2,I)=-FQ1(1,I)*ACY
       FQWK(3,I)= FQ1(2,I)
      ENDDO
      TMP1= FQ1(1,8)+FQ1(1,5)*RAB
      TMP2= FQ1(2,8)+FQ1(2,5)*RAB
       FQWK(1,9)=-TMP1*AQX
       FQWK(2,9)=-TMP1*ACY
       FQWK(3,9)= TMP2
C
      DO 110 I=1,4
         R01(1,I)= R01(1,I)+FQWK(1,1)*WORK(I+5)
         R01(2,I)= R01(2,I)+FQWK(2,1)*WORK(I+5)
         R01(3,I)= R01(3,I)+FQWK(3,1)*WORK(I+5)
  110 CONTINUE
      DO 120 I=5,8
         R01(1,I)= R01(1,I)+FQWK(1,2)*WORK(I+1)
         R01(2,I)= R01(2,I)+FQWK(2,2)*WORK(I+1)
         R01(3,I)= R01(3,I)+FQWK(3,2)*WORK(I+1)
  120 CONTINUE
      DO 130 I=9,12
         R01(1,I)= R01(1,I)+FQWK(1,3)*WORK(I-3)
         R01(2,I)= R01(2,I)+FQWK(2,3)*WORK(I-3)
         R01(3,I)= R01(3,I)+FQWK(3,3)*WORK(I-3)
  130 CONTINUE
      DO 140 I=13,16
         R01(1,I)= R01(1,I)+FQWK(1,4)*WORK(I-7)
         R01(2,I)= R01(2,I)+FQWK(2,4)*WORK(I-7)
         R01(3,I)= R01(3,I)+FQWK(3,4)*WORK(I-7)
  140 CONTINUE
      DO 150 I=17,20
         R01(1,I)= R01(1,I)+FQWK(1,5)*WORK(I-11)
         R01(2,I)= R01(2,I)+FQWK(2,5)*WORK(I-11)
         R01(3,I)= R01(3,I)+FQWK(3,5)*WORK(I-11)
  150 CONTINUE
      DO 160 I=21,25
         R01(1,I)= R01(1,I)+FQWK(1,6)*WORK(I-20)
         R01(2,I)= R01(2,I)+FQWK(2,6)*WORK(I-20)
         R01(3,I)= R01(3,I)+FQWK(3,6)*WORK(I-20)
  160 CONTINUE
      DO 170 I=26,30
         R01(1,I)= R01(1,I)+FQWK(1,7)*WORK(I-25)
         R01(2,I)= R01(2,I)+FQWK(2,7)*WORK(I-25)
         R01(3,I)= R01(3,I)+FQWK(3,7)*WORK(I-25)
  170 CONTINUE
      DO 180 I=31,35
         R01(1,I)= R01(1,I)+FQWK(1,8)*WORK(I-30)
         R01(2,I)= R01(2,I)+FQWK(2,8)*WORK(I-30)
         R01(3,I)= R01(3,I)+FQWK(3,8)*WORK(I-30)
  180 CONTINUE
      DO 190 I=36,40
         R01(1,I)= R01(1,I)+FQWK(1,9)*WORK(I-35)
         R01(2,I)= R01(2,I)+FQWK(2,9)*WORK(I-35)
         R01(3,I)= R01(3,I)+FQWK(3,9)*WORK(I-35)
  190 CONTINUE
C
      XMD4XY=XMD4*AQXY
      XMD4X=-XMD4*AQX
      XMD4Y=-XMD4*ACY
      DO 210 I=1,5
         R02(1,I)= R02(1,I)+(FQ2(1,I)*AQX2-FQ1(1,I))*XMD4
         R02(2,I)= R02(2,I)+(FQ2(1,I)*ACY2-FQ1(1,I))*XMD4
         R02(3,I)= R02(3,I)+(FQ2(3,I)     -FQ1(1,I))*XMD4
         R02(4,I)= R02(4,I)+FQ2(1,I)*XMD4XY
         R02(5,I)= R02(5,I)+FQ2(2,I)*XMD4X
         R02(6,I)= R02(6,I)+FQ2(2,I)*XMD4Y
  210 CONTINUE
C
      TMP1= FQ2(1,8)+FQ2(1,5)*RAB
      TMP2= FQ2(2,8)+FQ2(2,5)*RAB
      TMP3= FQ2(3,8)+FQ2(3,5)*RAB
      TMP4= FQ1(1,8)+FQ1(1,5)*RAB
      FQWK(1,5)= TMP1*AQX2-TMP4
      FQWK(2,5)= TMP1*ACY2-TMP4
      FQWK(3,5)= TMP3     -TMP4
      FQWK(4,5)= TMP1*AQXY
      FQWK(5,5)=-TMP2*AQX
      FQWK(6,5)=-TMP2*ACY
      DO I=6,9
      FQWK(1,I)= FQ2(1,I)*AQX2-FQ1(1,I)
      FQWK(2,I)= FQ2(1,I)*ACY2-FQ1(1,I)
      FQWK(3,I)= FQ2(3,I)     -FQ1(1,I)
      FQWK(4,I)= FQ2(1,I)*AQXY
      FQWK(5,I)=-FQ2(2,I)*AQX
      FQWK(6,I)=-FQ2(2,I)*ACY
      ENDDO
      DO 220 I=6,9
       DO 220 J=1,6
         R02(J,I)= R02(J,I)+FQWK(J,6)*WORK(I)
  220 CONTINUE
      DO 230 I=10,13
       DO 230 J=1,6
         R02(J,I)= R02(J,I)+FQWK(J,7)*WORK(I-4)
  230 CONTINUE
      DO 240 I=14,17
       DO 240 J=1,6
         R02(J,I)= R02(J,I)+FQWK(J,8)*WORK(I-8)
  240 CONTINUE
      DO 250 I=18,21
       DO 250 J=1,6
         R02(J,I)= R02(J,I)+FQWK(J,5)*WORK(I-12)
  250 CONTINUE
      DO 260 I=22,26
       DO 260 J=1,6
         R02(J,I)= R02(J,I)+FQWK(J,9)*WORK(I-21)
  260 CONTINUE
C
      FQ3(1,1)= FQ3(1,4)+FQ3(1,1)*RAB
      FQ3(2,1)= FQ3(2,4)+FQ3(2,1)*RAB
      FQ3(3,1)= FQ3(3,4)+FQ3(3,1)*RAB
      FQ3(4,1)= FQ3(4,4)+FQ3(4,1)*RAB
      FQ2(1,5)= FQ2(1,8)+FQ2(1,5)*RAB
      FQ2(2,5)= FQ2(2,8)+FQ2(2,5)*RAB
C
      DO 310 I=1,4
         R03( 1,I)= R03( 1,I)+(FQ3(1,I)*AQX2-FQ2(1,I+4)*THREE)*XMD4X
         R03( 2,I)= R03( 2,I)+(FQ3(1,I)*AQX2-FQ2(1,I+4)      )*XMD4Y
         R03( 3,I)= R03( 3,I)+(FQ3(2,I)*AQX2-FQ2(2,I+4)      )*XMD4
         R03( 4,I)= R03( 4,I)+(FQ3(1,I)*ACY2-FQ2(1,I+4)      )*XMD4X
         R03( 5,I)= R03( 5,I)+(FQ3(2,I)*AQXY                 )*XMD4
         R03( 6,I)= R03( 6,I)+(FQ3(3,I)     -FQ2(1,I+4)      )*XMD4X
         R03( 7,I)= R03( 7,I)+(FQ3(1,I)*ACY2-FQ2(1,I+4)*THREE)*XMD4Y
         R03( 8,I)= R03( 8,I)+(FQ3(2,I)*ACY2-FQ2(2,I+4)      )*XMD4
         R03( 9,I)= R03( 9,I)+(FQ3(3,I)     -FQ2(1,I+4)      )*XMD4Y
         R03(10,I)= R03(10,I)+(FQ3(4,I)     -FQ2(2,I+4)*THREE)*XMD4
  310 CONTINUE
      WORK2(1)=-(FQ3(1,5)*AQX2-FQ2(1,9)*THREE)*AQX
      WORK2(2)=-(FQ3(1,5)*AQX2-FQ2(1,9)      )*ACY
      WORK2(3)= (FQ3(2,5)*AQX2-FQ2(2,9)      )
      WORK2(4)=-(FQ3(1,5)*ACY2-FQ2(1,9)      )*AQX
      WORK2(5)= (FQ3(2,5)*AQXY               )
      WORK2(6)=-(FQ3(3,5)     -FQ2(1,9)      )*AQX
      WORK2(7)=-(FQ3(1,5)*ACY2-FQ2(1,9)*THREE)*ACY
      WORK2(8)= (FQ3(2,5)*ACY2-FQ2(2,9)      )
      WORK2(9)=-(FQ3(3,5)     -FQ2(1,9)      )*ACY
      WORK2(10)= (FQ3(4,5)     -FQ2(2,9)*THREE)
      DO 320 I=5,8
       DO 320 J=1,10
         R03( J,I)= R03( J,I)+WORK2(J)*WORK(I+1)
  320 CONTINUE
C
      AQX4= AQX2*AQX2
      ACY4= ACY2*ACY2
      X2Y2= AQX2*ACY2
      R04( 1)= R04( 1)+(FQ4(1)*AQX4-FQ3(1,5)*SIX*AQX2
     *                 +FQ2(1,9)*THREE)*XMD4
      R04( 2)= R04( 2)+(FQ4(1)*AQX2-FQ3(1,5)*THREE)*XMD4XY
      R04( 3)= R04( 3)+(FQ4(2)*AQX2-FQ3(2,5)*THREE)*XMD4X
      R04( 4)= R04( 4)+(FQ4(1)*X2Y2-FQ3(1,5)*AQX2-FQ3(1,5)*ACY2
     *                 +FQ2(1,9))*XMD4
      R04( 5)= R04( 5)+(FQ4(2)*AQX2-FQ3(2,5)      )*XMD4Y
      R04( 6)= R04( 6)+(FQ4(3)*AQX2-FQ3(1,5)*AQX2-FQ3(3,5)
     *                 +FQ2(1,9))*XMD4
      R04( 7)= R04( 7)+(FQ4(1)*ACY2-FQ3(1,5)*THREE)*XMD4XY
      R04( 8)= R04( 8)+(FQ4(2)*ACY2-FQ3(2,5)      )*XMD4X
      R04( 9)= R04( 9)+(FQ4(3)     -FQ3(1,5)      )*XMD4XY
      R04(10)= R04(10)+(FQ4(4)     -FQ3(2,5)*THREE)*XMD4X
      R04(11)= R04(11)+(FQ4(1)*ACY4-FQ3(1,5)*SIX*ACY2
     *                 +FQ2(1,9)*THREE)*XMD4
      R04(12)= R04(12)+(FQ4(2)*ACY2-FQ3(2,5)*THREE)*XMD4Y
      R04(13)= R04(13)+(FQ4(3)*ACY2-FQ3(1,5)*ACY2-FQ3(3,5)
     *                 +FQ2(1,9))*XMD4
      R04(14)= R04(14)+(FQ4(4)     -FQ3(2,5)*THREE)*XMD4Y
      R04(15)= R04(15)+(FQ4(5)     -FQ3(3,5)*SIX+FQ2(1,9)*THREE)*XMD4
      RETURN
      END
