C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, ROOT, FRGINF
C 23 JUL 04 - MWS - UNLINDEP: ORTHONORMALIZE LINEARLY INDEPENDENT SPACE
C 19 MAY 04 - DGF - SAVE RELATIVISTIC KINETIC ENERGY FOR DK (VIRIAL BUG)
C 01 APR 04 - DGF - CHANGES TO ADD THE FMO METHOD
C 16 JAN 04 - DGF - HFPVP,PVPVINT: FIX RESC LOOP BALANCING
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  4 NOV 03 - MWS - ONLINE IN GAMESS, REPLACING RESC.SRC
C  1 MAY 02 - DGF,TN - ADD NUMERIC DK GRADIENT
C  1 MAY 02 - DGF - PARALLELISE PVPINT, DPVPINT,HFPVP: FIX PARALLEL BUGS
C 23 APR 02 - DGF - ALLOW CHOOSING DOUBLE/QUADRUPLE PRECISION IN RIS
C 23 NOV 01 - TN  - ADD THE DOUGLAS-KROLL (DK) METHOD
C  6 SEP 01 - MWS - RESCX,RESCGX: ALLOW FOR CHECK RUNS
C 13 JUN 00 - DGF,TN - IMPLEMENT INT. UNCONTRACTED ENERGY AND GRADIENT
C 15 AUG 00 - DGF - PVPINT,DPVPINT,HFPVP: FIX L SHELL INTS, USE GLDIAG
C 25 MAR 00 - DGF - MAKE THE BASIS SET LINEAR DEPENDENCE CASE WORK
C  6 MAR 00 - DGF - CODE NESC GRADIENT
C 29 FEB 00 - DGF,TN - CODE RESC GRADIENT
C 10 JAN 00 - MWS - SYNCHRONIZE RELWFN COMMON BLOCK
C 21 DEC 99 - TN,DGF - ADD RESC METHOD (RELATIVISTIC CORRECTIONS)
C
C*MODULE QREL     *DECK PVPINT
      SUBROUTINE PVPINT(PVP2,LL2,DBUG,ISAVE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBUG,IANDJ,NORM,UNCON,ISEPS,USEPS,GOPARR,DSKWRK,MASWRK,NXT
      DOUBLE PRECISION MOROKM
C
      DIMENSION PVP2(LL2)
C
      PARAMETER (MXSH=5000,MXGTOT=20000,MXATM=2000,MAXDEN=25*MXATM,
     *           LENABC=2000,MXCHRM=1)
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /STV   / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
      COMMON /SYMIND/ TOL,II,JJ,LIT,LJT,MINI,MINJ,MAXI,MAXJ,IANDJ
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /ISEPS / ISEPS,USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION DIJ(225),PVP(225)
      DIMENSION IJX(35),IJY(35),IJZ(35)
      DIMENSION     XV(7,6,6),    YV(7,6,6),    ZV(7,6,6)
C     DIMENSION     XV(7,6,5),    YV(7,6,5),    ZV(7,6,5)
C                     (6,6,5) IS ACTUALLY NEEDED, BUT QRDDERIJ USES LD=7
      DIMENSION DXVDIJ(5,5,6),DYVDIJ(5,5,6),DZVDIJ(5,5,6)
C     DIMENSION DXVDIJ(5,5,5),DYVDIJ(5,5,5),DZVDIJ(5,5,5)
C
      DATA MAXRYS /9/
      DATA RLN10  /2.30258D+00/
      DATA ZERO   /0.0D+00/
      DATA ONE    /1.0D+00/
      DATA PI212  /1.1283791670955D+00/
      DATA SQRT3  /1.73205080756888D+00/
      DATA SQRT5  /2.23606797749979D+00/
      DATA SQRT7  /2.64575131106459D+00/
      DATA IJX    / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     1              4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     2              5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     3              3, 1, 3, 2, 2/
      DATA IJY    / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     1              1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     2              1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     3              1, 3, 2, 3, 2/
      DATA IJZ    / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     1              1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     2              1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     3              3, 3, 2, 2, 3/
      DATA MOROKM/8HMOROKUMA/,ANONE/8HNONE    /
C
      TOL=RLN10*ITOL
      NORM=NORMF.NE.1.OR.NORMP.NE.1
      UNCON=RMETHOD.NE.ANONE.AND.MOD(MODQR,2).EQ.1
      IF (RUNTYP.EQ.MOROKM) THEN
         CALL STINT1(ISTART,IEND,JSTART,LOCIJ,NATST,NATED,ISAVE,L1,L2)
      ELSE
         ISTART = 1
         IEND   = NSHELL
         JSTART = 1
         LOCIJ  = 0
         NATST  = 1
         NATED  = NAT+NCHMAT
         ISAVE  = 0
         L1 = NUM
         IF(UNCON) L1=NUMU
         L2 = (L1*(L1+1))/2
         IF(ISEPS) NATED = NAT+NPS
      END IF
C
      CALL VCLR(PVP2,1,LL2)
C
      NXT = IBTYP.EQ.1
      NEXT = -1
      MINE = -1
      ICOUNT=0
C
C     ----- CALCULATE -PVP- TERM -----
C
      IF(DBUG) WRITE(IW,9999)
C
C     ----- ISHELL -----
C
      DO 9000 II=ISTART,IEND
C
          IF(GOPARR.AND.NXT) THEN
            MINE = MINE + 1
            IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF(NEXT.NE.MINE) GO TO 9000
          END IF
C
      I=KATOM(II)
      XI=C(1,I)
      YI=C(2,I)
      ZI=C(3,I)
      I1=KSTART(II)
      I2=I1+KNG(II)-1
      LIT=KTYPE(II)
      MINI=KMIN(II)
      MAXI=KMAX(II)
      LOCI=KLOC(II)-MINI-LOCIJ
C
      LITDER=LIT+1
C
C     ----- JSHELL -----
C
      DO 8000 JJ=JSTART,II
        ICOUNT=ICOUNT+1
        IF(GOPARR.AND..NOT.NXT.AND.MOD(ICOUNT,NPROC).NE.ME) GOTO 8000
      J=KATOM(JJ)
      XJ=C(1,J)
      YJ=C(2,J)
      ZJ=C(3,J)
      J1=KSTART(JJ)
      J2=J1+KNG(JJ)-1
      LJT=KTYPE(JJ)
      MINJ=KMIN(JJ)
      MAXJ=KMAX(JJ)
      LOCJ=KLOC(JJ)-MINJ-LOCIJ
C
      LJTDER=LJT+1
C
      IANDJ=II.EQ.JJ
C     IANDJ=.FALSE.
      RR=(XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      NROOTS=(LITDER+LJTDER-2)/2+1
      IF(NROOTS.GT.MAXRYS) THEN
         WRITE(IW,9997) MAXRYS,LITDER,LJTDER,NROOTS
         STOP 1234
      ENDIF
C
      IJ=(MAXI-MINI+1)*(MAXJ-MINJ+1)
      CALL VCLR(PVP,1,IJ)
C
C     ----- I PRIMITIVE -----
C
      DO 7000 IG=I1,I2
      AI=EX(IG)
      ARRI=AI*RR
      AXI=AI*XI
      AYI=AI*YI
      AZI=AI*ZI
      CSI=CS(IG)
      CPI=CP(IG)
      CDI=CD(IG)
      CFI=CF(IG)
      CGI=CG(IG)
C
C     ----- J PRIMITIVE -----
C
      JGMAX=J2
      DO 6000 JG=J1,JGMAX
      AJ=EX(JG)
      AA=AI+AJ
      AA1=ONE/AA
      DUM=AJ*ARRI*AA1
      IF(DUM.GT.TOL) GO TO 6000
      FAC= EXP(-DUM)
      CSJ=CS(JG)
      CPJ=CP(JG)
      CDJ=CD(JG)
      CFJ=CF(JG)
      CGJ=CG(JG)
      AX=(AXI+AJ*XJ)*AA1
      AY=(AYI+AJ*YJ)*AA1
      AZ=(AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR -----
C
      IJ=0
      DO 360 I=MINI,MAXI
      GO TO (110,120,220,220,130,220,220,140,220,220,
     1       150,220,220,160,220,220,220,220,220,170,
     2       180,220,220,190,220,220,220,220,220,200,
     3       220,220,210,220,220),I
  110 DUM1=CSI*FAC
      GO TO 220
  120 DUM1=CPI*FAC
      GO TO 220
  130 DUM1=CDI*FAC
      GO TO 220
  140 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 220
  150 DUM1=CFI*FAC
      GO TO 220
  160 IF(NORM) DUM1=DUM1*SQRT5
      GO TO 220
  170 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 220
  180 DUM1=CGI*FAC
      GO TO 220
  190 IF(NORM) DUM1=DUM1*SQRT7
      GO TO 220
  200 IF(NORM) DUM1=DUM1*SQRT5/SQRT3
      GO TO 220
  210 IF(NORM) DUM1=DUM1*SQRT3
  220 CONTINUE
C
      JMAX=MAXJ
      IF(IANDJ) JMAX=I
      DO 360 J=MINJ,JMAX
      GO TO (230,250,350,350,260,350,350,270,350,350,
     1       280,350,350,290,350,350,350,350,350,300,
     2       310,350,350,320,350,350,350,350,350,330,
     3       350,350,340,350,350),J
  230 DUM2=DUM1*CSJ
      GO TO 350
  250 DUM2=DUM1*CPJ
      GO TO 350
  260 DUM2=DUM1*CDJ
      GO TO 350
  270 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 350
  280 DUM2=DUM1*CFJ
      GO TO 350
  290 IF(NORM) DUM2=DUM2*SQRT5
      GO TO 350
  300 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 350
  310 DUM2=DUM1*CGJ
      GO TO 350
  320 IF(NORM) DUM2=DUM2*SQRT7
      GO TO 350
  330 IF(NORM) DUM2=DUM2*SQRT5/SQRT3
      GO TO 350
  340 IF(NORM) DUM2=DUM2*SQRT3
  350 CONTINUE
C
      IJ=IJ+1
  360 DIJ(IJ)=DUM2
C
C     ----- PVP -----
C
      AAX=AA*AX
      AAY=AA*AY
      AAZ=AA*AZ
      DO 500 KAT=NATST,NATED
      ZNUC=-ZAN(KAT)
      CX=C(1,KAT)
      CY=C(2,KAT)
      CZ=C(3,KAT)
      XX=AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
      IF(NROOTS.LE.3) CALL RT123
      IF(NROOTS.EQ.4) CALL ROOT4
      IF(NROOTS.EQ.5) CALL ROOT5
      IF (NROOTS.GE.6) CALL ROOT6
      DO 420 IROOT=1,NROOTS
      UU=U(IROOT)*AA
      WW=W(IROOT)*ZNUC
      TT=ONE/(AA+UU)
      T = SQRT(TT)
      X0=(AAX+UU*CX)*TT
      Y0=(AAY+UU*CY)*TT
      Z0=(AAZ+UU*CZ)*TT
C
      DO 410 J=1,LJTDER
      NJ=J
      DO 410 I=1,LITDER
      NI=I
      CALL STVINT
      XV(I,J,IROOT)=XINT
      YV(I,J,IROOT)=YINT
      ZV(I,J,IROOT)=ZINT*WW
  410 CONTINUE
C
      CALL QRDDERIJ
     $           (DXVDIJ(1,1,IROOT),DYVDIJ(1,1,IROOT),DZVDIJ(1,1,IROOT),
     1                XV(1,1,IROOT),    YV(1,1,IROOT),    ZV(1,1,IROOT),
     2            LIT,LJT,AI,AJ,5)
C
  420 CONTINUE
C
      IJ=0
      DO 450 I=MINI,MAXI
      IX=IJX(I)
      IY=IJY(I)
      IZ=IJZ(I)
      JMAX=MAXJ
      IF(IANDJ) JMAX=I
      DO 440 J=MINJ,JMAX
      JX=IJX(J)
      JY=IJY(J)
      JZ=IJZ(J)
      IJ=IJ+1
C
C     DUMII STANDS FOR PI*V*PI INTEGRAL
C
      DUMXX=ZERO
      DUMYY=ZERO
      DUMZZ=ZERO
      DO 434 K=1,NROOTS
        DUMXX=DUMXX + DXVDIJ(IX,JX,K)* YV   (IY,JY,K)* ZV   (IZ,JZ,K)
        DUMYY=DUMYY +  XV   (IX,JX,K)*DYVDIJ(IY,JY,K)* ZV   (IZ,JZ,K)
        DUMZZ=DUMZZ +  XV   (IX,JX,K)* YV   (IY,JY,K)*DZVDIJ(IZ,JZ,K)
  434 CONTINUE
      DUM=(DUMXX+DUMYY+DUMZZ)*(AA1*PI212)
      PVP(IJ)=PVP(IJ)+DUM*DIJ(IJ)
  440 CONTINUE
  450 CONTINUE
C
  500 CONTINUE
C
 6000 CONTINUE
 7000 CONTINUE
C
C     ----- SET UP PVP MATRICES -----
C
      IJ=0
      DO 7500 I=MINI,MAXI
      JMAX=MAXJ
      IF(IANDJ) JMAX=I
      DO 7500 J=MINJ,JMAX
      IJ=IJ+1
      NN=((LOCI+I)*(LOCI+I-1))/2+(LOCJ+J)
      PVP2(NN)= PVP(IJ)
 7500 CONTINUE
 8000 CONTINUE
 9000 CONTINUE
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
      IF (GOPARR) CALL DDI_GSUMF(910,PVP2,LL2)
C
      RETURN
 9999 FORMAT(/,10X,20(1H-),/,10X,'1 ELECTRON PVP INTEGRALS',
     2       /,10X,20(1H-))
 9997 FORMAT(' IN -PVPINT- THE RYS QUADRATURE IS NOT IMPLEMENTED',
     1       ' BEYOND -NROOTS- = ',I2,/,
     2       ' LIT,LJT,NROOTS= ',3I3)
      END
C
C*MODULE QREL     *DECK QRDDERIJ
      SUBROUTINE QRDDERIJ(DDXDIJ,DDYDIJ,DDZDIJ,X,Y,Z,LIT,LJT,AI,AJ,LD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(7,*),Y(7,*),Z(7,*),DDXDIJ(LD,1),DDYDIJ(LD,1),
     *          DDZDIJ(LD,1)
C
C     ----- APPLY MOMENTUM P ON BOTH BRA AND KET (FOR PVP INTEGRALS)
C
      DDXDIJ(1,1)= X(2,2)*((AJ+AJ)*(AI+AI))
      DDYDIJ(1,1)= Y(2,2)*((AJ+AJ)*(AI+AI))
      DDZDIJ(1,1)= Z(2,2)*((AJ+AJ)*(AI+AI))
C
      IF(LJT.EQ.1) GO TO 330
      DO 320 J=2,LJT
      DDXDIJ(1,J)=(X(2,J+1)*(AJ+AJ)-X(2,J-1)*(J-1))*(AI+AI)
      DDYDIJ(1,J)=(Y(2,J+1)*(AJ+AJ)-Y(2,J-1)*(J-1))*(AI+AI)
      DDZDIJ(1,J)=(Z(2,J+1)*(AJ+AJ)-Z(2,J-1)*(J-1))*(AI+AI)
  320 CONTINUE
  330 IF(LIT.EQ.1) RETURN
      DO 340 I=2,LIT
      DDXDIJ(I,1)=(X(I+1,2)*(AI+AI)-X(I-1,2)*(I-1))*(AJ+AJ)
      DDYDIJ(I,1)=(Y(I+1,2)*(AI+AI)-Y(I-1,2)*(I-1))*(AJ+AJ)
      DDZDIJ(I,1)=(Z(I+1,2)*(AI+AI)-Z(I-1,2)*(I-1))*(AJ+AJ)
  340 CONTINUE
      IF(LJT.EQ.1) RETURN
      DO 350 I=2,LIT
        DO 350 J=2,LJT
          DDXDIJ(I,J)=(X(I+1,J+1)*(AJ+AJ)-X(I+1,J-1)*(J-1))*(AI+AI)
     2               -(X(I-1,J+1)*(AJ+AJ)-X(I-1,J-1)*(J-1))*(I-1)
          DDYDIJ(I,J)=(Y(I+1,J+1)*(AJ+AJ)-Y(I+1,J-1)*(J-1))*(AI+AI)
     2               -(Y(I-1,J+1)*(AJ+AJ)-Y(I-1,J-1)*(J-1))*(I-1)
          DDZDIJ(I,J)=(Z(I+1,J+1)*(AJ+AJ)-Z(I+1,J-1)*(J-1))*(AI+AI)
     2               -(Z(I-1,J+1)*(AJ+AJ)-Z(I-1,J-1)*(J-1))*(I-1)
  350 CONTINUE
      RETURN
      END
C
C*MODULE QREL     *DECK RESCX
      SUBROUTINE RESCX(H0,SAO,TAO,DWINT,L1,L2,DBUG,ISAVE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      DIMENSION H0(*),SAO(*),TAO(*),DWINT(*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /FMCOM / X(1)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      DATA TWO/2.0D+00/
      DATA CHECK/8HCHECK   /
C
C     ----- RELATIVISTIC SCHEME BY ELIMINATING SMALL-COMPONENTS -----
C                (USE POTENTIAL, NORMALIZED, SYMMETRIZED HAMILTONIAN)
C
C        CODED BY TAKAHITO NAKAJIMA AND PORTED BY DMITRI FEDOROV.
C        DEPARTMENT OF APPLIED CHEMISTRY, UNIVERSITY OF TOKYO (1998-99).
C
C     INTERESTINGLY, RESC DRIVER DOES NOT NEED TO KNOW IF IT IS RUNNING
C     IN THE INTERNALLY UNCONTRACTED MODE (SINCE THE CONTRACTION IS
C     PERFORMED ELSEWHERE (ONEEI)).
C
      CALL DERCHK(NDER)
      L2=(L1*(L1+1))/2
      L3=L1*L1
      CALL VALFM(LOADFM)
      LVEXT = LOADFM + 1
      LTAO2 = LVEXT + L2
      LXX = LTAO2 + L2
      LSD = LXX + L3
      LUU = LSD + L1
      LTD = LUU + L3
      LEI = LTD + L1
      LAI = LEI + L1
      LQI = LAI + L1
      LXU = LQI + L1
      LSXU = LXU + L3
      LYY =  LSXU + L3
      LXMM = LYY + L3
      LTWRK0 = LXMM + L3
      LTWRK1A = LTWRK0 + L1
      LTWRK2A = LTWRK1A + L2
      LTWRK2B = LTWRK2A + L3
      LGIVEB = LTWRK2B + L3
      LINDB = LGIVEB + L1*8
      LAST = LINDB + (L1 -1)/NWDVAR+1
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     ----- CALCULATE VEXT AND TAO*2 -----
C
      CALL VSUB(TAO,1,H0,1,X(LVEXT),1,L2)
      CALL DCOPY(L2,TAO,1,X(LTAO2),1)
      CALL DSCAL(L2,TWO,X(LTAO2),1)
      IF(NDER.GT.0) THEN
         CALL DAWRIT(IDAF,IODA,X(LTAO2),L2,NDARELB+4,0)
         CALL DAWRIT(IDAF,IODA,DWINT,L2,NDARELB+5,0)
         CALL DAWRIT(IDAF,IODA,X(LVEXT),L2,NDARELB+6,0)
         CALL DAWRIT(IDAF,IODA,SAO,L2,NDARELB+20,0)
C
C        WE MUST RECORD THE OVERLAPS SINCE THEY CAN BE IN THE
C        UNCONTRACTED BASIS AND RECORD 12 ALWAYS HAS THEM IN
C        THE CONTRACTED BASIS.
C
      ENDIF
C
C     RI STYLE:
C     0 - ORIGINAL RESC (ALL CARTESIAN FUNCTIONS IN RI)
C     1 - HONDO (REMOVE LINEAR DEPENDENT FUNCTIONS FROM RI)
C
      KAMERAD=MOD(MODQR/2,2)
C
C     ----- CALCULATE RELATIVISTIC HAMILTONIAN -----
C
      CALL ERESC(H0,SAO,X(LTAO2),X(LVEXT),DWINT,X(LXX),X(LSD),X(LUU),
     *           X(LTD),X(LEI),X(LXU),X(LSXU),X(LYY),X(LXMM),
     *           X(LTWRK0),X(LTWRK1A),X(LTWRK2A),X(LTWRK2B),
     *           X(LGIVEB),X(LINDB),NDER,KAMERAD,L0,L1,L2)
C
C     ----- WRITE RELATIVISTIC HAMILTONIAN TO DAF -----
C     H0 AT THIS POINT CAN BE IN THE UNCONTRACTED BASIS SET.
C     THE HAMILTONIAN WILL THEN HAVE TO BE REWRITTEN IN THE
C     CONTRACTED BASIS.
C
      CALL DAWRIT(IDAF,IODA,H0,L2,11,0)
      IF(ISAVE.EQ.1) CALL DAWRIT(IDAF,IODA,H0,L2,311,0)
C       FOR MOROKUMA DECOMPOSITION
      IF(DBUG) THEN
         WRITE(IW,*)'$$$$$ DARWIN INTEGRALS $$$$$'
         CALL PRTRI(DWINT,L1)
         WRITE(IW,*)'$$$$$ RELATIVISTIC H0 $$$$$'
         CALL PRTRIL(H0,L1)
      ENDIF
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C
C*MODULE QREL     *DECK QRSXU
      SUBROUTINE QRSXU(SXU,SAO,TAO2,XX,SD,UU,TD,XU,TWRK0,TWRK1A,
     $                 GIVEB,INDB,NDER,KAMERAD,L1,L0,L2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL UNCON,QUADRE,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /BASSPH/ QMTTOL,ISPHER
      DIMENSION SAO(*),TAO2(*),XX(L1,*),SD(*),UU(L1,*),TD(*),XU(L1,*),
     *          SXU(L1,*),TWRK0(*),TWRK1A(*),GIVEB(L1,*),INDB(*)
C
C     --- CONSTRUCT SXU WHICH WILL BE USED TO GET M = (SXU) Y (SXU)T
C
C     --- CONSTRUCT ORTHONORMAL(ON) BASIS FUNCTIONS (GET V,SD ABD X) ---
C
      CALL DCOPY(L2,SAO,1,TWRK1A,1)
C
      QMTSAV=QMTTOL
      QMTTOL=QRQMT
      UNCON =MOD(MODQR  ,2).EQ.1
      QUADRE=MOD(MODQR/8,2).EQ.1
C
C     MESSY POINT ABOUT ISPHER: UNCON=.FALSE. WITH KAMERAD.EQ.1 CALLS
C     QMTSYM, WHICH KEEPS SPHERICAL FUNCTIONS IF ISPHER.NE.1, SO WE
C     DO THE SAME FOR UNCON=.TRUE.
C
      IF(KAMERAD.EQ.0.OR.(KAMERAD.EQ.1.AND.ISPHER.NE.1.AND.UNCON)) THEN
C        ORIGINAL RESC PAPER'S METHODOLOGY
         CALL ONBAS(TWRK1A,XX,SD,GIVEB,INDB,TWRK0,NDER,L1,L0,QMTTOL,
     *              QUADRE)
      ELSE IF(KAMERAD.EQ.1) THEN
         L3=L1*L1
         IF(UNCON) THEN
C
C           WE NOW GENERATE PSEUDO-SALC-MATRIX. IT WILL HAVE SPHERICAL
C           CONTAMINANTS ELIMINATED BUT NO SYMMETRY BUILT-IN. THIS
C           MATRIX WILL BE USED FOR THE RIS. THE REASON FOR ALL THAT
C           IS THE SCHIZOPHRENIAC IDEA TO WORK WITH NORMAL AND
C           INTERNALLY UNCONTRACTED BASES AT THE SAME TIME. THIS
C           MEANS THAT SYMMETRY SHELL-RELATED INFORMATION IN COMMON
C           BLOCKS SUCH AS MAPCTR IS NOT STORED FOR BOTH SO WE CANNOT
C           GENERATE NORMAL Q-MATRIX. SXU, XX AND XU ARE USED AS
C           TEMPORARY STORAGE.
C
            CALL SPHBCKU(XU,L1,L1,NSALCU)
            IF(NPRINT.NE.-5.AND.MASWRK) WRITE(IW,9000) L1-NSALCU,NSALCU
C           TRANSFORM TO SPHERICAL HARMONICS
            CALL TFTRI(TWRK1A,SAO,XU,TWRK0,NSALCU,L1,L1)
C           PROCEED WITH ELIMINATION AS USUAL IN SPH. HARM. SPACE
            CALL ONBAS(TWRK1A,SXU,SD,GIVEB,INDB,TWRK0,NDER,NSALCU,L0,
     *                 QMTTOL,QUADRE)
C           BACK-TRANSFORM TO CARTESIANS
            CALL MRARBR(XU,L1,L1,NSALCU,SXU,NSALCU,L0,XX,L1)
            CALL TFTRI(TWRK1A,SAO,XX,TWRK0,L0,L1,L1)
C           WRITE(6,*) 'DOUBLE CHECK: SHOULD BE A UNIT MATRIX S'''
C           CALL PRTRI(TWRK1A,L0)
         ELSE
            CALL QMTSYM(SAO,UU,XX,SD,GIVEB,TWRK0,L0,L1,L2,L3,.FALSE.)
         ENDIF
      ENDIF
C
C     ----- TRANSFORM KINETIC ENERGY MATRIX -----
C
      CALL TFTRI(TWRK1A,TAO2,XX,TWRK0,L0,L1,L1)
C
C     ----- DIAGONALIZE KINETIC ENERGY MATRIX (GET U AND TD) -----
C        CAUTION: UU IS SAVED AS L1*L0 ARRAY EVEN THOUGH IT IS L0*L0
C
C--   CALL VCLR(UU,1,L1*L0)
C--   LENH = (L0*L0+L0)/2
C--   KORDER =0
C--   CALL EVVRSP(IW,L0,L0,LENH,L1,TWRK1A,GIVEB,INDB,TD,UU,KORDER,IERR)
      IERR=0
      IF(QUADRE) THEN
C        KINETIC ENERGY MATRIX HAS LARGE VALUES, SO THRESHOLD IS LARGE
C        CALL DCLEAN((L0*L0+L0)/2,TWRK1A,1,1.0D-12)
         CALL QKIN2(XX,UU,TD,L0,L1,L2)
C        CALL DQDIAG(TWRK1A,UU,TD,L1,L0,0)
      ELSE
         CALL GLDIAG(L1,L0,L0,TWRK1A,GIVEB,TD,UU,IERR,INDB)
      ENDIF
C
      IF(IERR.NE.0) THEN
         WRITE(IW,*) 'T DIAG FAILED WITH GLDIAG, TRYING JACOBI',IERR
C        CALL JACDG(TWRK1A,UU,TD,GIVEB(1,1),GIVEB(1,2),L1,L0)
         CALL ABRT
      ENDIF
C
      IF(NDER.GT.0) THEN
         CALL DAWRIT(IDAF,IODA,TD,L0,   NDARELB+13,0)
         CALL DAWRIT(IDAF,IODA,UU,L1*L0,NDARELB+14,0)
C        CAUTION: UU IS SAVED AS L1*L0 ARRAY EVEN THOUGH IT IS L0*L0
      ENDIF
C
C     ----- CALCULATE X*U -----
C     ----- CALCULATE S*XU -----
C
      CALL MRARBR(XX,L1,L1,L0,UU,L1,L0,XU,L1)
      CALL MTARBR(SAO,L1,XU,L0,SXU,L1,1)
      QMTTOL=QMTSAV
      RETURN
C
 9000 FORMAT(1X,'SPHERICAL CONTAMINANTS ELIMINATED:',I5,'; KEPT',I6,
     *       ' ORBITALS')
      END
C
C
C*MODULE QREL     *DECK QRMV
      SUBROUTINE QRMV(H0,XMM,YY,EI,SXU,TWRK1A,TWRK2A,L1,L0,L2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     ----- CALCULATE MASS-VELOCITY (MV) TERM -----
C
      DIMENSION H0(*),XMM(L1,*),YY(L1,*),EI(*),SXU(L1,*),TWRK1A(*),
     *          TWRK2A(L1,*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
C
C     -- CONSTRUCT RELATIVISTIC KINETIC ENERGY M IN COORDINATE SPACE --
C     M = (SXU) EI (SXU)T, EI IS KINETIC ENERGY IN P-SPACE, A.K.A. EP
C
      CALL VCLR(YY,1,L1*L1)
      DO I=1,L0
C       SHIFT -MC**2
        YY(I,I)=EI(I)-CLIG2
      ENDDO
C
C     ------- CALCULATE (SXU)EI*(SXU)T
C
      CALL MRARBR(SXU,L1,L1,L0,YY,L1,L0,TWRK2A,L1)
      CALL MRARTR(TWRK2A,L1,L1,L0,SXU,L1,L1,XMM,L1)
C     WRITE(6,*) 'HAVING HQR',YY(1,1),XMM(1,1),SXU(1,1)
C
      CALL CPYSQT(XMM,TWRK1A,L1,1)
      CALL VADD(H0,1,TWRK1A,1,H0,1,L2)
C
C     ----- WRITE MV TERM TO DAF -----
C
      CALL DAWRIT(IDAF,IODA,H0,L2,13,0)
C     WRITE(6,*) 'QRMV: SAVING REL KIN',TWRK1A(1),TWRK1A(2)
C
      RETURN
      END
C
C*MODULE QREL     *DECK ONBAS
      SUBROUTINE ONBAS(SS,V,E,GIVEB,INDB,TWRK0,NDER,L1,L0,TOL,QUADRE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK,QUADRE
      DIMENSION SS(*),V(L1,*),E(*),GIVEB(L1,*),INDB(*),TWRK0(*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      DATA ZERO/0.0D+00/,ONE/1.0D+00/,TEN9/1.0D+09/,SMALL/1.0D-11/
C
C     ----- GENERATE A SET OF ORTHONORMAL ORBITALS WITH
C           NEAR LINEAR-DEPENDENCIES REMOVED.
C
C--   CALL VCLR(V,1,L1*L1)
C--   IERR=0
C--   CALL JACDG(SS,V,E,GIVEB(1,1),GIVEB(1,2),L1,L1)
C
C--   CALL VCLR(V,1,L1*L1)
C--   CALL GIVEIS(L1,L1,L1,SS,GIVEB,INDB,E,V,IERR)
      IERR=0
      IF(QUADRE) THEN
C            OVERLAP HAS SMALL (<=1) VALUES, SO THE THRESHOLD IS SMALL
         CALL DCLEAN((L1*L1+L1)/2,SS,1,1.0D-14)
         CALL DQDIAG(SS,V,E,L1,L1,0)
      ELSE
         CALL GLDIAG(L1,L1,L1,SS,GIVEB,E,V,IERR,INDB)
      ENDIF
      IF(IERR.NE.0) THEN
         WRITE(IW,*) 'S DIAG FAILED',IERR
         CALL ABRT
      ENDIF
      CALL DCOPY(L1,E,1,TWRK0,1)
C
C     ----- ELIMINATE EIGENVECTORS FOR WHICH EIGENVALUE IS LESS -----
C         THAN -TOL- ( DEFINED IN /QMTOPT/ AND DEFAULT = 1.0E-04)
C     THE LINEARLY DEPENDENT COMPONENTS ARE MOVED TO THE BACK. THEY ARE
C     NEEDED FOR THE GRADIENT AND ARE NOT USED FOR THE ENERGY.
C     THIS ALSO DESTROYS THE SORTED ORDER OF E. DOES IT MATTER?
C
      ESMALL=TEN9
      NSMALL=0
      ISALC =1
      JSALC =L1
  520 CONTINUE
        IF(E(ISALC).LT.ESMALL) ESMALL=E(ISALC)
        E(ISALC)=ONE/ SQRT(E(ISALC))
        IF(ISALC.GE.JSALC) GOTO 530
        IF(TWRK0(ISALC).GE.TOL) THEN
          ISALC=ISALC+1
        ELSE
          NSMALL=NSMALL+1
          CALL DSWAP(L1,V(1,ISALC),1,V(1,JSALC),1)
          CALL DSWAP( 1,E(ISALC),  1,E(JSALC),1)
          CALL DSWAP( 1,TWRK0(ISALC),1,TWRK0(JSALC),1)
          JSALC=JSALC-1
        ENDIF
      GOTO 520
  530 CONTINUE
C
      L0=L1-NSMALL
      IF(MASWRK.AND.  .NOT.(NPRINT.EQ.-5.OR. NPRINT.EQ.-2183)
     *      .AND.  NSMALL.NE.0)
     *   WRITE(IW,9999) ESMALL,NSMALL,TOL,L0
      NQMTR=L0
      IF(NDER.GT.0) THEN
C        IF(NSMALL.NE.0) THEN
C           WRITE(IW,9000)
C           CALL ABRT
C        ENDIF
         CALL DAWRIT(IDAF,IODA,TWRK0,L1,NDARELB+7,0)
         CALL DAWRIT(IDAF,IODA,V, L1*L1,NDARELB+8,0)
      ENDIF
C
C     ----- FORM CANONICAL ORTHONORMAL ORBITALS -----
C
      DO J=1,L1
         EJ = E(J)
         DO I=1,L1
            QIJ=V(I,J)
            QIJ=QIJ*EJ
            IF( ABS(QIJ).LT.SMALL) QIJ=ZERO
            V(I,J)=QIJ
         ENDDO
      ENDDO
      IF(NDER.GT.0) CALL DAWRIT(IDAF,IODA,V,L1*L1,NDARELB+12,0)
C
      RETURN
 9999 FORMAT(
     1 ' THE SMALLEST EIGENVALUE OF THE OVERLAP MATRIX IS ',1P,E10.2/
     2 ' THERE ARE ',I5,' OVERLAP MATRIX EIGENVALUES LESS THAN ',
     *      1P,E10.2/
     3 ' THE NUMBER OF LINEARLY INDEPENDENT ORBITALS KEPT IS',I6/1X)
      END
C
C*MODULE QREL     *DECK QREI
      SUBROUTINE QREI(EI,TD,L0)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION EI(*),TD(*)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
C
C     ----- CALCULATE EI, KINETIC ENERGY IN MOMENTUM SPACE-----
C
      CLIG4=CLIG2*CLIG2
      DO I=1,L0
        EI(I)= SQRT(CLIG2*TD(I)+CLIG4)
      ENDDO
C
      RETURN
      END
C
C*MODULE QREL     *DECK ERESC
      SUBROUTINE ERESC(H0,SAO,TAO2,VEXT,DWINT,XX,SD,UU,TD,EI,XU,SXU,YY,
     $                 XMM,TWRK0,TWRK1A,TWRK2A,TWRK2B,GIVEB,INDB,
     $                 NDER,KAMERAD,L0,L1,L2)
C
C     ----- RELATIVISTIC SCHEME TO ELIMINATE THE SMALL-COMPONENT -----
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H0(*),SAO(*),TAO2(*),VEXT(*),DWINT(*),XX(L1,*),
     *          SD(*),UU(L1,*),TD(*),EI(*),XU(L1,*),SXU(L1,*),
     *          YY(L1,*),XMM(L1,*),TWRK0(*),TWRK1A(*),
     *          TWRK2A(L1,*),TWRK2B(L1,*),GIVEB(L1,8),INDB(L1)
C
C     ----- CONSTRUCT SXU -----
C     THIS CALL WILL ALSO SET THE VALUE OF -L0-
C
      CALL QRSXU(SXU,SAO,TAO2,XX,SD,UU,TD,XU,TWRK0,TWRK1A,GIVEB,INDB,
     $           NDER,KAMERAD,L1,L0,L2)
C
      CALL VCLR(H0,1,L2)
C
C     ------- CALCULATE EI
C
      CALL QREI(EI,TD,L0)
C
C     ------- CALCULATE SQR TERM
C
      CALL QRMV(H0,XMM,YY,EI,SXU,TWRK1A,TWRK2A,L1,L0,L2)
C
C     ADD RESC TRANSFORMED PVP CONTRIBUTION
      CALL QRESC(H0,XMM,YY,TD,SXU,XU,DWINT,TWRK1A,TWRK2A,TWRK2B,
     *           L1,L0,L2,1)
C     ADD RESC TRANSFORMED VEXT CONTRIBUTION
      CALL QRESC(H0,XMM,YY,TD,SXU,XU,VEXT, TWRK1A,TWRK2A,TWRK2B,
     *           L1,L0,L2,2)
C     FOR KINETIC ONLY DEBUG
C     CALL VADD(H0,1,VEXT,1,H0,1,L2)
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRESC
      SUBROUTINE QRESC(H0,XMM,YY,TD,SXU,XU,RELINT,TWRK1A,TWRK2A,
     $                 TWRK2B,L1,L0,L2,ITYPE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H0(*),XMM(L1,*),YY(L1,*),SXU(L1,*),XU(L1,*),RELINT(*),
     *          TD(*),TWRK1A(*),TWRK2A(L1,*),TWRK2B(L1,*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
C
C     ----- CALCULATE DARWIN-LIKE TERM -----
C
C        (SXU) * YY-LEFT * (XU)T
C
      CALL QRYY(YY,YY,TD,L1,L0,.TRUE.,.FALSE.,ITYPE)
      CALL MRARBR(SXU,L1,L1,L0,YY,L1,L0,TWRK2A,L1)
      CALL MRARTR(TWRK2A,L1,L1,L0,XU,L1,L1,XMM,L1)
C
C     FOR SPIN-ORBIT (ITYPE==1) AND GRADIENT SAVE THE TRANSFORMATION
C
      NDARE=NDARELB+(ITYPE-1)*2
      CALL DAWRIT(IDAF,IODA,XMM,L1*L1,NDARE,0)
C
C     ----- MULTIPLY ONE-ELECTRON INTEGRAL V -----
C
      CALL CPYTSQ(RELINT,TWRK2A,L1,1)
      CALL MRARBR(XMM,L1,L1,L1,TWRK2A,L1,L1,TWRK2B,L1)
C
C        (SXU) * YY-RIGHT * (XU)T
C
      CALL QRYY(YY,YY,TD,L1,L0,.FALSE.,.TRUE.,ITYPE)
      CALL MRARBR(SXU,L1,L1,L0,YY,L1,L0,TWRK2A,L1)
      CALL MRARTR(TWRK2A,L1,L1,L0,XU,L1,L1,XMM,L1)
C
C     FOR SPIN-ORBIT (ITYPE==1) AND GRADIENT SAVE THE TRANSFORMATION
C
      CALL DAWRIT(IDAF,IODA,XMM,L1*L1,NDARE+1,0)
C
C     ----- (SXU) * YY-LEFT * (XU)T * W * (XU) * YY-RIGHT * (SXU)T
C           WHERE W IS EITHER PVP (ITYPE==1) OR V (ITYPE==2).
C
      CALL MRARTR(TWRK2B,L1,L1,L1,XMM,L1,L1,TWRK2A,L1)
C
C     ----- SYMMETRIZED HAMILTONIAN -----
C
      CALL SYMTRZE(TWRK2A,L1,L1)
C
C     ----- CONSTRUCT H0 -----
C
      CALL CPYSQT(TWRK2A,TWRK1A,L1,1)
C     FOR KINETIC ONLY DEBUG COMMENT OUT THE NEXT LINE
      CALL VADD(H0,1,TWRK1A,1,H0,1,L2)
C
      RETURN
      END
C*MODULE QREL     *DECK RESCGX
      SUBROUTINE RESCGX(DENAB,L1,L2,L3)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW,NXT,SVGPAR,UNCON
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      DIMENSION DENAB(L2),LDVREC(2),LPOT(2),LXMML(2),LXMMR(2)
      DATA CHECK/8HCHECK   /
C
C     __________ RESC ENERGY GRADIENT _____________
C
C     CODED BY DMITRI FEDOROV AND TAKAHITO NAKAJIMA.
C
C     DEPARTMENT OF APPLIED CHEMISTRY, UNIVERSITY OF TOKYO (1999-2000).
C
      UNCON=MOD(MODQR,2).EQ.1
      CALL VALFM(LOADFM)
      LV= LOADFM + 1
      LPVP=LV+L2
      LSD= LPVP+L2
      LTD= LSD+L1
      LVV= LTD+L1
      LXX= LVV+L3
      LUU= LXX+L3
      LTT= LUU+L3
      LSS= LTT+L2
      LPP= LSS+L2
      LSX= LPP+L3
      LXU= LSX+L3
      LSXU=LXU+L3
      LDSD=LSXU+L3
      LDTD=LDSD+L1
      LDSINT=LDTD+L1
      LDTINT=LDSINT+L2
      LDVINT=LDTINT+L2
      LDX=LDVINT+L2
      LDU=LDX+L3
      LWORK=LDU+L3
      LWRK1A=LWORK+L1
      LWRK2A=LWRK1A+L2
      LWRK2B=LWRK2A+L3
      LWRK2C=LWRK2B+L3
      LWRK2D=LWRK2C+L3
      LYYDL=LWRK2D+L3
      LXMML(1)=LYYDL+L3
      LXMMR(1)=LXMML(1)+L3
      LXMML(2)=LXMMR(1)+L3
      LXMMR(2)=LXMML(2)+L3
      LDYYR=LXMMR(2)+L3
      LYYR=LDYYR+L3
      LYYDR=LYYR+L3
      LDPP=LYYDR+L3
      LDYYL=LDPP+L3
      LYYL=LDYYL+L3
      LTMPL=LYYL+L3
      LTMPR=LTMPL+L3
      LAST= LTMPR+L3
C     NOTE THAT GRADEN REQUIRES SPECIAL STRUCTURE OF THE ARGUMENTS.
C     ALSO THERE IS SOME OVERLAP WITH PREVIOUS STORAGE
      LUNCDER=LWRK1A
      IF(UNCON) THEN
        LCONDER=LAST
        LU=LCONDER+(NUM*NUM+NUM)/2
        LAST=LU+L1*NUM
      ELSE
        LCONDER=LUNCDER
        LU=LAST
      ENDIF
C
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     S   OVERLAPS
C     YL1 LEFT TRANSFORMATION FOR PVP: SXU * YL * (XU)T
C     YR1 RIGHT TRANSFORMATION FOR PVP:SXU * YR * (XU)T
C     YL2 LEFT TRANSFORMATION FOR V:   SXU * YL * (XU)T
C     YR2 RIGHT TRANSFORMATION FOR V:  SXU * YR * (XU)T
C     T   NON-RELATIVISTIC KINETIC ENERGY *2
C     PVP RELATIVISTIC CORRECTION TO POTENTIAL ENERGY
C     V   NON-RELATIVISTIC POTENTIAL ENERGY
C     SD  EIGENVALUES OF OVERLAP
C     V   EIGENVECTORS OF OVERLAP
C     X   ORTHOGONALISATION TRANSFORMATIOB OF AOS
C     TD  EIGENVALUES OF DOUBLED KINETIC ENERGY
C     U   EIGENVECTORS OF KINETIC ENERGY
C     UB  TRANSFORMATION TO CONTRACTED BASIS SET
C
C     CALL DAREAD(IDAF,IODA,X(LSS),L2,12,0)
      CALL DAREAD(IDAF,IODA,X(LSS),L2,NDARELB+20,0)
      L0=NQMTR
C
      CALL DAREAD(IDAF,IODA,X(LXMML(1)),L3,NDARELB,0)
      CALL DAREAD(IDAF,IODA,X(LXMMR(1)),L3,NDARELB+1,0)
      CALL DAREAD(IDAF,IODA,X(LXMML(2)),L3,NDARELB+2,0)
      CALL DAREAD(IDAF,IODA,X(LXMMR(2)),L3,NDARELB+3,0)
      CALL DAREAD(IDAF,IODA,X(LTT),L2,   NDARELB+4,0)
      CALL DAREAD(IDAF,IODA,X(LPVP),L2,  NDARELB+5,0)
      CALL DAREAD(IDAF,IODA,X(LV),L2,    NDARELB+6,0)
      CALL DAREAD(IDAF,IODA,X(LSD),L1,   NDARELB+7,0)
      CALL DAREAD(IDAF,IODA,X(LVV),L1*L1,NDARELB+8,0)
      CALL DAREAD(IDAF,IODA,X(LXX),L1*L1,NDARELB+12,0)
      CALL DAREAD(IDAF,IODA,X(LTD),L0,   NDARELB+13,0)
      CALL DAREAD(IDAF,IODA,X(LUU),L1*L0,NDARELB+14,0)
      IF(UNCON) CALL DAREAD(IDAF,IODA,X(LU),L1*NUM,NDARELB+19,0)
C
C     ----- CALCULATE SX, XU, SXU -----
C
      CALL MTARBR(X(LSS),L1,X(LXX),L0,X(LSX),L1,1)
      CALL MRARBR(X(LXX),L1,L1,L0,X(LUU),L1,L0,X(LXU),L1)
      CALL MTARBR(X(LSS),L1,X(LXU),L0,X(LSXU),L1,1)
C
      SVDSKW=DSKWRK
      DSKWRK=.TRUE.
C
C     INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      NEXT = -1
      MINE = -1
C
      LRECS=0
      LRECT=3*NAT
      LRECPVP=2*3*NAT
      LRECV=3*3*NAT
      LREC=0
      DO IATM=1,NAT
        DO IXYZ=1,3
          LREC=LREC+1
          IF(GOPARR.AND.NXT) THEN
            MINE = MINE + 1
            IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF(NEXT.NE.MINE) GO TO 500
          END IF
          IF(GOPARR.AND..NOT.NXT.AND.MOD(LREC,NPROC).NE.ME) GOTO 500
C
C         READ DERIVATIVE OF ORBITAL OVERLAPS AND KINETIC ENERGY
C         NB: THROUGHOUT D/DA MEANS NUCLEAR COORDINATE "A" DERIVATIVE
C
          CALL RAREAD(NQRDAF,X(MORDA),X(LDSINT),L2,LRECS+LREC,0)
          CALL RAREAD(NQRDAF,X(MORDA),X(LDTINT),L2,LRECT+LREC,0)
C
C         DISGRACEFUL PART:
          SVGPAR=GOPARR
          GOPARR=.FALSE.
C
C         THE INTEGRAL TRANSFORMATIONS ARE DIVIDED BETWEEN NODES AS A
C         WHOLE, SO THAT EACH NODE DOES ITS OWN TRANSFORMATION
C         SEQUENTIALLY ALL CALLS TO TFTRI AND SUCH MUST BE UNGOPARRED
C         (INSANE AS IT IS).
C         NOTE THAT ALL READS FROM FILES MUST BE WITH PRISTINE GOPARR
C         QRDV AND QRDX OPERATE IN FULL SPACE (L1 VECTORS) AND QRDU IN
C         L0 SPACE.  DX/DA IS COMPUTED AS L1XL1 AND USED AS L1XL0.
C         THIS IS POSSIBLE SINCE LINERALY DEPENDENT VECTORS WERE PUSHED
C         TOWARD THE END.
C         QRDU HAS NO PROBLEM AS IT WORKS IN L0XL0 SPACE SO THAT NO
C         RECTANGULAR MATRIX INVERSION PROBLEM OCCURS.
C
C     ----- CALCULATE VT DV/DA AND DSD/DA -----
C
          CALL QRDV(X(LVV),X(LSD),X(LDSD),X(LWRK1A),X(LDSINT),X(LWRK2C),
     *              X(LWORK),L1)
C
C     ----- CALCULATE DX/DA -----
C
          CALL QRDX(X(LVV),X(LSD),X(LDSD),X(LWRK2C),X(LWRK2B),X(LWRK2A),
     *              X(LDX),L1)
C
C     ----- CALCULATE DU/DA -----
C
          CALL QRDU(X(LUU),X(LTD),X(LDTD),X(LTT),X(LXU),X(LDX),X(LWRK1A)
     $             ,X(LDTINT),X(LWRK2A),X(LWRK2B),X(LWRK2C),X(LWRK2D),
     *              X(LWORK),X(LDU),L1,L0)
C
          GOPARR=SVGPAR
C
C     ----- FIRST DERIVATIVE OF RELATIVISTIC KINETIC ENERGY -----
C
          CALL QRDMV(DE(IXYZ,IATM),DENAB,X(LSX),X(LXU),X(LSXU),X(LSS),
     *               X(LUU),X(LDTD),X(LPP),X(LDPP),X(LTD),X(LWORK),
     *               X(LDX),X(LDU),X(LDSINT),X(LUNCDER),X(LCONDER),
     *               X(LU),X(LWORK),X(LWRK2A),X(LWRK2B),X(LWRK2C),
     *               L1,L0,UNCON)
C
C     ----- FIRST DERIVATIVE OF RELATIVISTIC POTENTIAL ENERGY -----
C
C     ITYPE=1   ! PVP
C     ITYPE=2   ! VEXT (THAT IS, V=1/R)
C
          LPOT(1)=LPVP
          LPOT(2)=LV
          LDVREC(1)=LREC+LRECPVP
          LDVREC(2)=LREC+LRECV
C     FOR KINETIC ONLY DEBUG CHANGE TO 2,2
          DO ITYPE=1,2
C
            CALL CPYTSQ(X(LPOT(ITYPE)),X(LPP),L1,1)
C
C           READ DERIVATIVE PVP OR V INTEGRALS
C
            CALL RAREAD(NQRDAF,X(MORDA),X(LDVINT),L2,LDVREC(ITYPE),0)
C
C           COMPUTE DERIVATIVE OF RELATIVISTIC POTENTIAL ENERGY
C
            ITYPEX = ITYPE
            CALL QRDPOT(DE(IXYZ,IATM),DENAB,X(LPP),X(LXU),X(LSS),X(LXX),
     *                  X(LUU),X(LDTD),X(LDPP),X(LYYL),X(LYYR),X(LYYDL),
     *                  X(LYYDR),X(LTD),X(LDX),X(LDU),X(LXMML(ITYPEX)),
     *                  X(LXMMR(ITYPEX)),X(LDSINT),X(LDVINT),X(LWRK2A),
     *                  X(LDYYL),X(LWRK2C),X(LTMPL),X(LTMPR),X(LDYYR),
     *                  X(LUNCDER),X(LCONDER),X(LU),X(LWORK),X(LWRK2B),
     *                  L1,L0,UNCON,ITYPEX)
          ENDDO
  500   CONTINUE
        ENDDO
      ENDDO
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
      DSKWRK=SVDSKW
C
C     WE MUST NOT GLOBAL SUM THE GRADIENT SINCE THIS IS DONE ELSEWHERE
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C
C*MODULE QREL     *DECK QRDV
      SUBROUTINE QRDV(VV,SD,DSD,TWRK1A,DSINT,VTDV,WORK,L1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION SD(*),VV(L1,*),DSD(*),TWRK1A(*),DSINT(*),VTDV(L1,*),
     *          WORK(*)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DATA ZERO /0.0D+00/,TINY/1.0D-13/
C
C     DERIVATIVE OF V
C
      SMALL=QRTOL
C
C     ----- VT DS/DA V -----
C
      CALL TFTRI(TWRK1A,DSINT,VV,WORK,L1,L1,L1)
C
C     ----- DSD/DA -----
C
      DO I=1,L1
         II=(I*(I+1))/2
         DSD(I)=TWRK1A(II)
      ENDDO
C
C     ----- VT DV/DA -----
C
      KL=0
      DO K=1,L1
        DO L=1,K
          KL=KL+1
          IF(K.NE.L) THEN
            TEMP=SD(L)-SD(K)
C           IF(ABS(TEMP).GE.SMALL) THEN
            IF(ABS(TEMP/MAX(SD(L),SD(K))).GE.SMALL.AND.
     *         ABS(TEMP).GT.TINY) THEN
              VTDV(K,L)= TWRK1A(KL)/TEMP
              VTDV(L,K)=-TWRK1A(KL)/TEMP
            ELSE
              VTDV(K,L)=ZERO
              VTDV(L,K)=ZERO
C             IF(ABS(TWRK1A(KL)).LT.SMALL) THEN
C               VTDV(K,L)=1.0D+00
C               VTDV(L,K)=1.0D+00
C             ENDIF
            ENDIF
          ELSE
            VTDV(K,K)=ZERO
          ENDIF
        ENDDO
      ENDDO
C
      RETURN
      END
C
C
C*MODULE QREL     *DECK QRDX
      SUBROUTINE QRDX(VV,SD,DSD,VTDV,TWRK2B,TWRK2A,DX,L1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VV(L1,*),SD(*),DSD(*),VTDV(L1,*),TWRK2A(L1,*),
     *          TWRK2B(L1,*),DX(L1,*)
      DATA ONE/1.0D+00/,HALF/0.5D+00/
C
C     DERIVATIVE OF X = V*SQRT(SD)
C
C     ------- DV/DA = V VT DV/DA -------
C
      CALL MRARBR(VV,L1,L1,L1,VTDV,L1,L1,TWRK2B,L1)
C
C     ------- DV/DA S**-1/2 -------
C
      CALL VCLR(TWRK2A,1,L1*L1)
      DO I=1,L1
        TWRK2A(I,I)=ONE/ SQRT(SD(I))
      ENDDO
      CALL MRARBR(TWRK2B,L1,L1,L1,TWRK2A,L1,L1,DX,L1)
C
C     ----- -1/2 V S**-3/2 DSD/DA -----
C
      CALL VCLR(TWRK2A,1,L1*L1)
      DO I=1,L1
        TWRK2A(I,I)=ONE/(SD(I)* SQRT(SD(I)))*DSD(I)
      ENDDO
      CALL MRARBR(VV,L1,L1,L1,TWRK2A,L1,L1,TWRK2B,L1)
      CALL DAXPY(L1*L1,-HALF,TWRK2B,1,DX,1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRDU
      SUBROUTINE QRDU(UU,TD,DTD,TT,XU,DX,TWRK1A,DTINT,TWRK2A,TWRK2B,
     *                TWRK2C,TWRK2D,WORK,DU,L1,L0)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DIMENSION TD(*),DTD(*),UU(L1,*),TT(*),XU(L1,*),DX(L1,*),TWRK1A(*),
     *          DTINT(*),TWRK2D(L1,*),TWRK2A(L0,*),TWRK2B(L1,*),
     *          TWRK2C(L0,*),WORK(*),DU(L1,*)
      DATA ZERO /0.0D+00/,TWO/2.0D+00/,TINY/1.0D-14/
C
C     DERIVATIVE OF U
C
      SMALL=QRTOL
C
C     NOTE THE L0 DIMENSIONS ON TWRK2A AND TWRK2C
C
C     ----- (XU)T DT/DA (XU) -----
C     DOUBLING IS BECAUSE DTINT ARE DT/DA AND HERE WE NEED DERIVATIVE
C     OF P**2 = 2*T
C
      CALL TFTRI(TWRK1A,DTINT,XU,WORK,L0,L1,L1)
      CALL DSCAL((L0*(L0+1))/2,TWO,TWRK1A,1)
C
      CALL CPYTSQ(TWRK1A,TWRK2C,L0,1)
C
C     ----- DX/DA U -----
C
      CALL MRARBR(DX,L1,L1,L0,UU,L1,L0,TWRK2B,L1)
C
C     ----- (XU)T TAO2 (DX/DA*U) + (DX/DA*U)T TAO2 (XU) -----
C
      CALL MTARBR(TT,L1,TWRK2B,L0,TWRK2D,L1,1)
      CALL MRTRBR(XU,L1,L1,L0,TWRK2D,L1,L0,TWRK2A,L0)
      CALL ADDCC(L0,TWRK2A,L0,TWRK2C,L0,0)
C
C     ----- DTD/DA -----
C
      DO I=1,L0
         DTD(I)=TWRK2C(I,I)
      ENDDO
C
C     ----- UT DU/DA -----
C
      CALL VCLR(TWRK2A,1,L0*L0)
      DO K=1,L0
        DO L=1,L0
          IF(K.NE.L) THEN
            TEMP=TD(L)-TD(K)
C
C           IT IS A FINE POINT WHAT TO DO IN CASE OF NEARLY DEGENERATE
C           EIGENVALUES. THE DIFFERENCE BETWEEN S AND T IS IN THAT T
C           EIGENVALUES CAN BE VERY LARGE AND THUS IT IS IMPORTANT TO
C           LOOK FOR RELATIVE ERROR DUE TO THE LIMITED NUMBER OF
C           SIGNIFICANT FIGURES. S EIGENVALUES CAN BE VERY SMALL DUE
C           TO LINEAR DEGENERACIES. STILL USAGE OF RELATIVE ERROR SEEMS
C           BETTER??
C
C           IF(ABS(TEMP).GE.SMALL) THEN
            IF(ABS(TEMP/MAX(TD(L),TD(K))).GE.SMALL.AND.
     *         ABS(TEMP).GT.TINY) THEN
              TWRK2A(K,L)=TWRK2C(K,L)/TEMP
            ELSE
              TWRK2A(K,L)=ZERO
C             IF(ABS(TWRK2C(K,L)).LT.SMALL) TWRK2A(K,L)=1.0D+00
            ENDIF
          ELSE
            TWRK2A(K,K)=ZERO
          ENDIF
        ENDDO
      ENDDO
C
C     ----- DU/DA = U UT DU/DA -----
C
      CALL MRARBR(UU,L1,L0,L0,TWRK2A,L0,L0,DU,L1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRDMV
      SUBROUTINE QRDMV(EGRAD,DAB,SX,XU,SXU,SS,UU,DTD,XMV,XDMV,TD,EI,DX,
     *                 DU,DSINT,UNCDER,CONDER,UBAS,WORK,TWRK2A,TWRK2B,
     *                 TWRK2C,L1,L0,UNCON)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL UNCON
      DIMENSION DAB(*),SX(L1,*),XU(L1,*),SXU(L1,*),SS(*),UU(L1,*),DTD(*)
     *         ,XMV(L1,*),XDMV(L1,*),TD(*),EI(*),DX(L1,*),DU(L1,*),
     *          DSINT(*),UNCDER(*),CONDER(*),UBAS(L1,*),WORK(*),
     *          TWRK2A(L1,*),TWRK2B(L1,*),TWRK2C(L1,*)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DATA HALF /0.5D+00/
C
C     DERIVATIVE OF RELATIVISTIC KINETIC ENERGY
C
C     GIVEN BY (SXU) T (SXU)T, SECOND T STANDS FOR TRANSPOSED
C
      CALL VCLR(XMV,1,L1*L0)
      CALL QREI(EI,TD,L0)
      DO I=1,L0
        XMV(I,I)=EI(I)-CLIG2
      ENDDO
C
C     ----- SXU DT/DA (SXU)T -----
C
      CALL VCLR(TWRK2B,1,L1*L1)
      DO I=1,L0
        TWRK2B(I,I)=DTD(I)*HALF*CLIG2/EI(I)
      ENDDO
C     ALAS, POOR TRSQU. CALL TRSQU(TWRK2A,TWRK2B,SXU,TWRK0,L1,L0)
      CALL MRARBR(SXU,L1,L1,L0,TWRK2B,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,SXU,L1,L1,XDMV,L1)
C
C     ----- (DS/DA)XU T (SXU)T + C.C. -----
C
      CALL MTARBR(DSINT,L1,XU,L0,TWRK2B,L1,1)
      CALL MRARBR(TWRK2B,L1,L1,L0,XMV,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,SXU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,XDMV,L1,0)
C
C     ----- SXDU/DA T (SXU)T + C.C. -----
C
      CALL MRARBR(SX,L1,L1,L0,DU,L1,L0,TWRK2B,L1)
      CALL MRARBR(TWRK2B,L1,L1,L0,XMV,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,SXU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,XDMV,L1,0)
C
C     ----- SDX/DAU T (SXU)T + C.C. -----
C
      CALL MTARBR(SS,L1,DX,L0,TWRK2C,L1,1)
      CALL MRARBR(TWRK2C,L1,L1,L0,UU,L1,L0,TWRK2B,L1)
      CALL MRARBR(TWRK2B,L1,L1,L0,XMV,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,SXU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,XDMV,L1,0)
C
C     ----- MULTIPLY DENSITY MATRIX -----
C
      CALL CPYSQT(XDMV,UNCDER,L1,1)
C     EGRAD0=EGRAD
      CALL GRADEN(EGRAD,UNCDER,CONDER,UBAS,DAB,WORK,UNCON,1)
C     WRITE(6,*) 'GRAD DUE TO QRDMV',EGRAD-EGRAD0
      RETURN
      END
C
C*MODULE QREL     *DECK QRDPOT
      SUBROUTINE QRDPOT(EGRAD,DAB,PP,XU,SS,XX,UU,DTD,XDMV,YYL,YYR,YYDL,
     $                  YYDR,TD,DX,DU,XMML,XMMR,DSINT,DVINT,TWRK2A,DYYL,
     $                  TWRK2C,TMPL,TMPR,DYYR,UNCDER,CONDER,UBAS,WORK,
     *                  TWRK2B,L1,L0,UNCON,ITYPE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL UNCON
      DIMENSION DAB(*),PP(L1,*),XU(L1,*),SS(*),XX(L1,*),UU(L1,*),DTD(*),
     *          XDMV(L1,*),YYL(L1,*),YYR(L1,*),YYDL(L1,*),YYDR(L1,*),
     *          TD(*),DX(L1,*),DU(L1,*),XMML(L1,*),XMMR(L1,*),DSINT(*),
     *          DVINT(*),TWRK2A(L1,*),DYYL(L1,*),TWRK2C(L1,*),TMPL(L1,*)
     *         ,TMPR(L1,*),DYYR(L1,*),UNCDER(*),CONDER(*),UBAS(L1,*),
     *          WORK(*),TWRK2B(L1,*)
      DATA ONE/1.0D+00/
C
C     DERIVATIVE OF RELATIVISTIC POTENTIAL ENERGY E GIVEN BY
C
C     E = (H + (H)T)/2
C     H = SXU * YL * (XU)T * W * XU * YR * (SXU)T
C
C     W IS EITHER PVP (ITYPE=1) OR V (AKA VEXT) (ITYPE=2)
C
C     S(L1,L1) - OVERLAPS
C     X(L1,L0) - TRANSFORMATION TO ORTHOGONAL BASIS SET,
C                GIVEN BY (X)T*S*X=1
C     U(L0,L0) - EIGENVECTORS OF (DOUBLED) KINETIC ENERGY, P**2
C     YL(L0,L0) OQ  (ITYPE=1), O*SQRT(Q) (ITYPE=2), DIAGONAL MATRIX
C     YR(L0,L0) Q/O (ITYPE=1), SQRT(Q)/O (ITYPE=2), DIAGONAL MATRIX
C     FOR ITYPE=2 THE FACTOR OF 2MC IS PRESENTLY INCLUDED INTO YL
C
C     GRADIENT IS A DERIVATIVE OF A PRODUCT OF 13 MATRICES AND IT IS
C     GROUPED INTO 3 DERIVATIVE GRPS OF V (1), YL (6) AND YR (6) TERMS.
C
C     NOTE THAT THE SYMMETRISATION IS PERFORMED AT THE VERY END AND ALL
C     INTERMEDIATE C.C. COME FROM PRISTINE SYMMETRY IN THE ENERGY
C     EXPRESSION RESC GRADIENT IS *THE EXACT* ANALYTIC DERIVATIVE OF
C     RESC ENERGY
C
C     ----- LEFT AND RIGHT RESC TRANSFORMATION MATRICES (IN P-SPACE),
C                     YL AND YR AND THEIR DERIVATIVES
C
      CALL QRYY(YYL,YYR,TD,L1,L0,.TRUE.,.TRUE.,ITYPE)
      CALL QRDYY(DYYL,DYYR,TD,DTD,L1,L0,ITYPE)
C
C     DE/DA, PART I:  SXU * YL * (XU)T * DV/DA * XU * YR * (SXU)T
C
      CALL CPYTSQ(DVINT,TWRK2A,L1,1)
      CALL MRARBR(XMML,L1,L1,L1,TWRK2A,L1,L1,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L1,XMMR,L1,L1,XDMV,L1)
C
      CALL VCLR(TMPL,1,L1*L1)
      CALL VCLR(TMPR,1,L1*L1)
C
C     ------- X DU/DA YL (XU)T + C.C. (HERMITIAN CONJUGATE)
C
      CALL MRARBR(XX,L1,L1,L0,DU,L1,L0,TWRK2B,L1)
      CALL MRARBR(TWRK2B,L1,L1,L0,YYL,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,TMPL,L1,1)
C
C     ------- X DU/DA YR (XU)T + C.C.
C
      CALL MRARBR(TWRK2B,L1,L1,L0,YYR,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,TMPR,L1,1)
C
C     ------- DX/DA U YL (XU)T + C.C.
C
      CALL MRARBR(DX,L1,L1,L0,UU,L1,L0,TWRK2B,L1)
      CALL MRARBR(TWRK2B,L1,L1,L0,YYL,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,TMPL,L1,1)
C
C     ------- DX/DA U YR (XU)T + C.C.
C
      CALL MRARBR(TWRK2B,L1,L1,L0,YYR,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL ADDCC(L1,TWRK2A,L1,TMPR,L1,1)
C
C     ------- XU DYL/DA (XU)T
C
      CALL MRARBR(XU,L1,L1,L0,DYYL,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL DAXPY(L1*L1,ONE,TWRK2A,1,TMPL,1)
C
C     ------- XU DYR/DA (XU)T
C
      CALL MRARBR(XU,L1,L1,L0,DYYR,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL DAXPY(L1*L1,ONE,TWRK2A,1,TMPR,1)
C
C     S*{ (XDU/DAYL(XU)T + C.C.) + (DX/DAU*YL(XU)T + C.C.)
C         + XU DYL/DA (XU)T }
C     S*{ (XDU/DAYR(XU)T + C.C.) + (DX/DAU*YR(XU)T + C.C.)
C         + XU DYR/DA (XU)T }
C
      CALL MTARBR(SS,L1,TMPL,L1,YYDL,L1,1)
      CALL MTARBR(SS,L1,TMPR,L1,YYDR,L1,1)
C
C     ------- DS/DA XU YL (XU)T
C
      CALL MRARBR(XU,L1,L1,L0,YYL,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL MTARBR(DSINT,L1,TWRK2A,L1,TWRK2B,L1,1)
      CALL DAXPY(L1*L1,ONE,TWRK2B,1,YYDL,1)
C
C     ------- DS/DA XU YR (XU)T
C
      CALL MRARBR(XU,L1,L1,L0,YYR,L1,L0,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      CALL MTARBR(DSINT,L1,TWRK2A,L1,TWRK2B,L1,1)
      CALL DAXPY(L1*L1,ONE,TWRK2B,1,YYDR,1)
C
C     ----- YYDL V (XMMR)T -----
C     DE/DA, PART II: [ S*{ (XDU/DAYL(XU)T + C.C.)
C         + (DX/DAU*YL(XU)T + C.C.) + XU DYL/DA (XU)T }
C         + DS/DA XU YL (XU)T ] * V * [ SXU*YR*(XU)T ]T
C
      CALL MRARBR(YYDL,L1,L1,L1,PP,L1,L1,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L1,XMMR,L1,L1,TWRK2A,L1)
      CALL DAXPY(L1*L1,ONE,TWRK2A,1,XDMV,1)
C
C     ----- XMML V (YYDR)T -----
C     DE/DA, PART III: SXU * YL * (XU)T * V *
C       [ S*{ (X DU/DA YR(XU)T + C.C.) + (DX/DAU*YR(XU)T + C.C.)
C           + XU DYR/DA (XU)T } + DS/DA XU YR (XU)T ]T
C
      CALL MRARBR(XMML,L1,L1,L1,PP,L1,L1,TWRK2C,L1)
      CALL MRARTR(TWRK2C,L1,L1,L1,YYDR,L1,L1,TWRK2A,L1)
      CALL DAXPY(L1*L1,ONE,TWRK2A,1,XDMV,1)
C
C     FOR KINETIC ONLY DEBUG
C     CALL CPYTSQ(DVINT,XDMV,L1,1)
C
C     ----- MULTIPLY DENSITY MATRIX -----
C     (SYMMETRISE THE GRADIENT MATRIX FIRST)
C
C     EGRAD0=EGRAD
      CALL SYMSQT(L1,XDMV,L1,UNCDER,1)
      CALL GRADEN(EGRAD,UNCDER,CONDER,UBAS,DAB,WORK,UNCON,1)
C     WRITE(6,*) 'GRAD DUE TO QRDPOT',EGRAD-EGRAD0
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRYY
      SUBROUTINE QRYY(YYL,YYR,TD,L1,L0,LEFT,RIGHT,ITYPE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL LEFT,RIGHT
      DIMENSION YYL(L1,*),YYR(L1,*),TD(*)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      EQUIVALENCE (CLIG,C),(CLIG2,C2)
      DATA ONE  /1.0D+00/,TWO/2.0D+00/
C
C     COMPUTE THE LEFT AND RIGHT TRANSFORMATION MATRICES
C
C     YYL=O*Q           AND YYR=Q/O       (PVP, ITYPE=1)
C     YYL=2MC*SQRT(Q)*O AND YYR=SQRT(Q)/O (V,   ITYPE=2)
C
      IF(LEFT)  CALL VCLR(YYL,1,L1*L1)
      IF(RIGHT) CALL VCLR(YYR,1,L1*L1)
C
C     2MC DIVISION MUST BE CONSISTENT WITH QRDYY
C
      A=TWO*C
      B=ONE
C
      DO I=1,L0
        Q=ONE/(SQRT(TD(I)+C2)+C)
        COQ=SQRT(ONE+TD(I)*Q*Q)
        IF(ITYPE.EQ.1) THEN
           IF(LEFT)  YYL(I,I)=Q*Q*COQ/C
           IF(RIGHT) YYR(I,I)=C/COQ
        ELSE
           IF(LEFT)  YYL(I,I)=A*Q*SQRT(Q)*COQ/C
           IF(RIGHT) YYR(I,I)=B*C/COQ/SQRT(Q)
        ENDIF
      ENDDO
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRDYY
      SUBROUTINE QRDYY(YYL,YYR,TD,DTD,L1,L0,ITYPE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DIMENSION YYL(L1,*),YYR(L1,*),TD(*),DTD(*)
C
      EQUIVALENCE (CLIG,C),(CLIG2,C2)
C
      DATA ONE/1.0D+00/,TWO/2.0D+00/,THREE/3.0D+00/,FIVE/5.0D+00/
C
C     COMPUTE DERIVATIVES OF THE LEFT AND RIGHT TRANSFORMATION MATRICES
C
      CALL VCLR(YYL,1,L1*L1)
      CALL VCLR(YYR,1,L1*L1)
C
C     2MC DIVISION, MUST BE CONSISTENT WITH QRYY
C
      A=TWO*C
      B=ONE
C
      DO I=1,L0
        EP=SQRT(C2*TD(I)+C2*C2)
        Q=ONE/(SQRT(TD(I)+C2)+C)
        DQ=-Q*Q*C*DTD(I)/(TWO*EP)
        QDTD=Q*Q*DTD(I)
        COQ=SQRT(ONE+TD(I)*Q*Q)
        COQ3=COQ*COQ*COQ
        IF(ITYPE.EQ.1) THEN
         YYL(I,I)=((TWO+THREE*TD(I)*Q*Q)*DQ+Q*QDTD/TWO)*Q/C/COQ
         YYR(I,I)=-C*(TD(I)*Q*DQ+QDTD/TWO)/COQ3
        ELSE
         YYL(I,I)=A*((THREE+FIVE*TD(I)*Q*Q)*DQ+Q*QDTD)*SQRT(Q)/C/TWO/COQ
         YYR(I,I)=-B*C*((ONE/Q+THREE*TD(I)*Q)*DQ+QDTD)/TWO/COQ3/SQRT(Q)
        ENDIF
      ENDDO
C
      RETURN
      END
C*MODULE QREL     *DECK ADDCC
      SUBROUTINE ADDCC(N,A,LDA,B,LDB,MODE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,N),B(LDB,N)
C
C     ADD SQUARE MATRIX A AND ITS HERMITIAN CONJUGATE TO SQUARE MATRIX B
C     B = B + A + A-TRANSPOSE
C     MODE=0: LOWER TRIANGLE OF B IS USED AND THE RESULT COPIED TO UPPER
C     MODE=1: B IS USED AS IS
C
      DO I=1,N
        IF(MODE.EQ.0) THEN
          DO J=1,I
            B(I,J)=B(I,J)+A(I,J)+A(J,I)
            B(J,I)=B(I,J)
          ENDDO
        ELSE
          DO J=1,N
            B(I,J)=B(I,J)+A(I,J)+A(J,I)
          ENDDO
        ENDIF
      ENDDO
      RETURN
      END
C
C*MODULE QREL     *DECK DPVPINT
      SUBROUTINE DPVPINT(DENAB,DPVP,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DBG,NORM
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,SVDSKW
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      DIMENSION DIJ(225),IJX(35),IJY(35),IJZ(35),DPVP(L2,3,NAT),
     *          XV(7,6,6),YV(7,6,6),ZV(7,6,6),DXV(5,5,6),DYV(5,5,6),
     *          DZV(5,5,6),XVDIJ(6,5,6),YVDIJ(6,5,6),ZVDIJ(6,5,6),
     *          DXVDIJ(5,5,6),DYVDIJ(5,5,6),DZVDIJ(5,5,6),DENAB(*)
C    *          XV(7,6,5),YV(7,6,5),ZV(7,6,5),DXV(5,5,5),DYV(5,5,5),
C    *          DZV(5,5,5),XVDIJ(6,5,5),YVDIJ(6,5,5),ZVDIJ(6,5,5),
C    *          DXVDIJ(5,5,5),DYVDIJ(5,5,5),DZVDIJ(5,5,5),DENAB(*)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /FMCOM / X(1)
      COMMON /GRAD  / DE(3,MXATM)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,TWO=2.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00,SQRT7=2.64575131106459D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HTVDER   /,RESC/8HRESC    /
C
      DATA IJX / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *           4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *           5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *           3, 1, 3, 2, 2/
      DATA IJY / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *           1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *           1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *           1, 3, 2, 3, 2/
      DATA IJZ / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *           1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *           1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *           3, 3, 2, 2, 3/
C
C     ----- BASIS FUNCTION DERIVATIVE CONTRIBUTIONS TO GRADIENT -----
C                   INTEGRALS ARE OF TYPE <II'/PVP/JJ>
C          (THAT IS, < (DMU/DR) | PVP | NU> + DERIVATIVE OF KET)
C                NOTE THAT DPVP ARRAY IS PASSED FROM HFPVP
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      IF(DBG) WRITE(IW,9000)
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NDER=1
      FSC=ONE/CLIG
      FSC2=FSC*FSC/TWO
C
C INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = 0
      NEXT = -1
      MINE = -1
C
C     -D*ELMM- CONTAINS CONTRIBUTION TO THE FORCES ON THE
C     EXTERNAL CHARGES WHICH ARE DUE TO QM MOLECULE.
C
C     IF(NCHMAT.NE.0) THEN
C        DO 100 ICHARM=1,NCHMAT
C           DXELMM(ICHARM)=ZERO
C           DYELMM(ICHARM)=ZERO
C           DZELMM(ICHARM)=ZERO
C 100    CONTINUE
C     END IF
C
C     FIRST, LOOP OVER THE GENUINE QM MOLECULE (ICHARM=0 PASS),
C     AND THEN ANY EXTERNAL CHARGES WHICH CHARMM MAY BE USING.
C     THE LATTER MAKE A CONTRIBUTION TO THE GRADIENT OF THE ATOMS
C     IN THE QM MOLECULE, AS WELL AS TO THEMSELVES.
C
C     DO 2000 ICHARM=0,NCHMAT
C
C     ----- I SHELL
C
      DO 1400 II = 1,NSHELL
C
C           GO PARALLEL!
C
        IF(NXT .AND. GOPARR) THEN
           MINE = MINE + 1
           IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
           IF(NEXT.NE.MINE) GO TO 1400
        END IF
C
        IAT = KATOM(II)
        XI = C(1,IAT)
        YI = C(2,IAT)
        ZI = C(3,IAT)
        I1 = KSTART(II)
        I2 = I1+KNG(II)-1
        LIT = KTYPE(II)
        MINI = KMIN(II)
        MAXI = KMAX(II)
        LOCI = KLOC(II)-MINI
        LITDER = LIT + NDER + 1
C
C     ----- J SHELL
C
        DO 1300 JJ = 1,NSHELL
C
C           GO PARALLEL!
C
          IF((.NOT.NXT) .AND. GOPARR) THEN
             IPCOUNT = IPCOUNT + 1
             IF(MOD(IPCOUNT,NPROC).NE.ME) GO TO 1300
          END IF
C
          JAT = KATOM(JJ)
          XJ = C(1,JAT)
          YJ = C(2,JAT)
          ZJ = C(3,JAT)
          J1 = KSTART(JJ)
          J2 = J1+KNG(JJ)-1
          LJT = KTYPE(JJ)
          MINJ = KMIN(JJ)
          MAXJ = KMAX(JJ)
          LOCJ = KLOC(JJ)-MINJ
          LJTDER = LJT+1
          NROOTS = (LITDER+LJTDER-2)/2 + 1
          RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
C
C     ----- I PRIMITIVE
C
        DO 1100 IG = I1,I2
          AI = EX(IG)
          ARRI = AI*RR
          AXI = AI*XI
          AYI = AI*YI
          AZI = AI*ZI
          CSI=CS(IG)
          CPI=CP(IG)
          CDI=CD(IG)
          CFI=CF(IG)
          CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
          DO 1000 JG = J1,J2
            AJ = EX(JG)
            AA = AI+AJ
            AA1 = ONE/AA
            DUM = AJ*ARRI*AA1
            IF(DUM .GT. TOL) GO TO 1000
            FAC = EXP(-DUM)
            CSJ = CS(JG)
            CPJ = CP(JG)
            CDJ = CD(JG)
            CFJ = CF(JG)
            CGJ = CG(JG)
            AX = (AXI+AJ*XJ)*AA1
            AY = (AYI+AJ*YJ)*AA1
            AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
            IJ = 0
            DUM1 = ZERO
            DUM2 = DUM1
            DO 600 I=MINI,MAXI
              IF(I.EQ.1) DUM1=CSI*FAC
              IF(I.EQ.2) DUM1=CPI*FAC
              IF(I.EQ.5) DUM1=CDI*FAC
              IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.11) DUM1=CFI*FAC
              IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
              IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.21) DUM1=CGI*FAC
              IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
              IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
              IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
C
              DO 600 J = MINJ,MAXJ
                IF(J.EQ.1) DUM2=DUM1*CSJ
                IF(J.EQ.2) DUM2=DUM1*CPJ
                IF(J.EQ.5) DUM2=DUM1*CDJ
                IF(J.EQ.8.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.11) DUM2=DUM1*CFJ
                IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
                IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.21) DUM2=DUM1*CGJ
                IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
                IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
                IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
C
                IJ=IJ+1
                IF(RMETHOD.EQ.RESC) THEN
                  DEN=ONE
                  IF(LOCI+I.EQ.LOCJ+J) DEN=TWO
                ELSE
                  NN = IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                  DEN = FSC2/TWO*DENAB(NN)
C                 THE ABOVE "TWO" COMES FROM NESC ITSELF
                  DEN = DEN + DEN
                ENDIF
                DIJ(IJ)=DUM2*DEN
  600         CONTINUE
C
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
      DO 790 IC = 1,NAT
         ZNUC = -ZAN(IC)
         CX = C(1,IC)
         CY = C(2,IC)
         CZ = C(3,IC)
         XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
         IF(NROOTS.LE.3) CALL RT123
         IF(NROOTS.EQ.4) CALL ROOT4
         IF(NROOTS.EQ.5) CALL ROOT5
         IF(NROOTS.GE.6) CALL ROOT6
         DO 740 K = 1,NROOTS
            UU = AA*U(K)
            WW = W(K)*ZNUC
            TT = ONE/(AA+UU)
            T = SQRT(TT)
            X0 = (AAX+UU*CX)*TT
            Y0 = (AAY+UU*CY)*TT
            Z0 = (AAZ+UU*CZ)*TT
            DO 730 J = 1,LJTDER
               NJ = J
               DO 720 I = 1,LITDER
                  NI = I
                  CALL VINT
                  XV(I,J,K) = XINT
                  YV(I,J,K) = YINT
                  ZV(I,J,K) = ZINT*WW
  720          CONTINUE
  730       CONTINUE
C
C           GET <MU|PVP|NU>
C
            CALL QRDDERIJ(XVDIJ(1,1,K),YVDIJ(1,1,K),ZVDIJ(1,1,K),
     1         XV(1,1,K),YV(1,1,K),ZV(1,1,K),LIT+NDER,LJT,AI,AJ,5+NDER)
C
C           GET <(DMU/DR)|V|NU> AND <(DMU/DR)|PVP|NU>
C
            CALL GDERI(DXV(1,1,K),DYV(1,1,K),DZV(1,1,K),
     *                 XV(1,1,K),YV(1,1,K),ZV(1,1,K),7,LIT,LJT,AI)
            CALL DERI(DXVDIJ(1,1,K),DYVDIJ(1,1,K),DZVDIJ(1,1,K),
     *                XVDIJ(1,1,K),YVDIJ(1,1,K),ZVDIJ(1,1,K),LIT,LJT,AI)
C
  740    CONTINUE
         IJ=0
         DO 780 I=MINI,MAXI
           IX=IJX(I)
           IY=IJY(I)
           IZ=IJZ(I)
           DO 780 J=MINJ,MAXJ
             JX=IJX(J)
             JY=IJY(J)
             JZ=IJZ(J)
C     DUMIJI CONTAINS CONTRIBUTION TO (D/DRJ)*PI*V*PI IN LOOSE NOTATION
             DUMXXX=ZERO
             DUMXYX=ZERO
             DUMXZX=ZERO
             DUMYXY=ZERO
             DUMYYY=ZERO
             DUMYZY=ZERO
             DUMZXZ=ZERO
             DUMZYZ=ZERO
             DUMZZZ=ZERO
             DO 430 K=1,NROOTS
C                            <D/DRX PX V PX> <V> <V>
               DUMXXX=DUMXXX+DXVDIJ(IX,JX,K)* YV(IY,JY,K)* ZV(IZ,JZ,K)
C                            <PX V PX> <D/DRY V> <V>
               DUMXYX=DUMXYX+ XVDIJ(IX,JX,K)*DYV(IY,JY,K)* ZV(IZ,JZ,K)
C                            <PX V PX> <V> <D/DRZ V> ET CETERA
               DUMXZX=DUMXZX+ XVDIJ(IX,JX,K)* YV(IY,JY,K)*DZV(IZ,JZ,K)
               DUMYXY=DUMYXY+DXV(IX,JX,K)* YVDIJ(IY,JY,K)* ZV(IZ,JZ,K)
               DUMYYY=DUMYYY+ XV(IX,JX,K)*DYVDIJ(IY,JY,K)* ZV(IZ,JZ,K)
               DUMYZY=DUMYZY+ XV(IX,JX,K)* YVDIJ(IY,JY,K)*DZV(IZ,JZ,K)
               DUMZXZ=DUMZXZ+DXV(IX,JX,K)* YV(IY,JY,K)* ZVDIJ(IZ,JZ,K)
               DUMZYZ=DUMZYZ+ XV(IX,JX,K)*DYV(IY,JY,K)* ZVDIJ(IZ,JZ,K)
               DUMZZZ=DUMZZZ+ XV(IX,JX,K)* YV(IY,JY,K)*DZVDIJ(IZ,JZ,K)
  430        CONTINUE
             IJ=IJ+1
C
C         I.E. <MU|(D/DRJ)*P*V*P|NU>=SUM OVER I <MU|(D/DRJ)*PI*V*PI|NU>
C
             DUMX=DUMXXX+DUMYXY+DUMZXZ
             DUMY=DUMXYX+DUMYYY+DUMZYZ
             DUMZ=DUMXZX+DUMYZY+DUMZZZ
             DUMINT=DIJ(IJ)*AA1*PI212
             IF(RMETHOD.EQ.RESC) THEN
                NN = IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                DPVP(NN,1,IAT)=DPVP(NN,1,IAT)+DUMX*DUMINT
                DPVP(NN,2,IAT)=DPVP(NN,2,IAT)+DUMY*DUMINT
                DPVP(NN,3,IAT)=DPVP(NN,3,IAT)+DUMZ*DUMINT
             ELSE
                DE(1,IAT)=DE(1,IAT)+DUMX*DUMINT
                DE(2,IAT)=DE(2,IAT)+DUMY*DUMINT
                DE(3,IAT)=DE(3,IAT)+DUMZ*DUMINT
             ENDIF
  780    CONTINUE
  790 CONTINUE
C
 1000 CONTINUE
 1100 CONTINUE
 1300 CONTINUE
 1400 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1500,DPVP,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
        LREC=2*3*NAT
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),DPVP(1,J,I),L2,LREC,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,42(1H-)/10X,
     *      'GRADIENT INCLUDING PVP AO DERIVATIVE TERMS',/10X,42(1H-))
      END
C*MODULE QREL     *DECK HFPVP
      SUBROUTINE HFPVP(DENAB,DVINTS,L2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,DBG,NORM,GOPARR,DSKWRK,MASWRK,NXT
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      DIMENSION XV(7,6,6),YV(7,6,6),ZV(7,6,6),DENAB(*),DVINTS(L2,3,NAT),
     *          DXV(7,6,6),DYV(7,6,6),DZV(7,6,6),
     *          XVDIJ(5,5,6),YVDIJ(5,5,6),ZVDIJ(5,5,6),
     *          DXVDIJ(5,5,6),DYVDIJ(5,5,6),DZVDIJ(5,5,6),
     *          IJX(35),IJY(35),IJZ(35),DIJ(225)
C     DIMENSION XV(7,6,5),YV(7,6,5),ZV(7,6,5),DENAB(*),DVINTS(L2,3,NAT),
C    *          DXV(7,6,5),DYV(7,6,5),DZV(7,6,5),
C    *          XVDIJ(5,5,5),YVDIJ(5,5,5),ZVDIJ(5,5,5),
C    *          DXVDIJ(5,5,5),DYVDIJ(5,5,5),DZVDIJ(5,5,5),
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /GRAD  / DE(3,MXATM)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,TWO=2.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00)
      PARAMETER (SQRT7=2.64575131106459D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HHELFEY  /,RESC/8HRESC    /
      DATA IJX/ 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *          4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *          5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *          3, 1, 3, 2, 2/
      DATA IJY/ 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *          1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *          1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *          1, 3, 2, 3, 2/
      DATA IJZ/ 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *          1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *          1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *          3, 3, 2, 2, 3/
C
C     ----- CALCULATE HELLMANN-FEYNMAN PVP INTEGRALS -----
C               (THAT IS, <MU| P(DV/DR)P |NU>)
C
C     RESC NEEDS PURE INTEGRALS, NESC CONTRACTS THEM WITH DENSITY
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      IF(DBG) WRITE(IW,9000)
C
      IF(RMETHOD.EQ.RESC) CALL VCLR(DVINTS,1,L2*3*NAT)
      FSC=ONE/CLIG
      FSC2=FSC*FSC/TWO
      NDER=1
C     1 FOR V AND, IN ADDITION, 1 FOR BRA, 1 FOR KET = 3RD DERIVATIVE
C
      TOL = RLN10*ITOL
      NORM = NORMF.NE.1 .OR. NORMP.NE.1
C
C     INITIALIZATION FOR PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = 0
      NEXT = -1
      MINE = -1
C
C     ----- I SHELL
C
      DO 600 II = 1,NSHELL
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         MINE = MINE + 1
         IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.MINE) GO TO 600
      END IF
C
      I = KATOM(II)
      XI = C(1,I)
      YI = C(2,I)
      ZI = C(3,I)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      LITDER=LIT+1
C
C     ----- J SHELL
C
      DO 580 JJ = 1,II
C
C           GO PARALLEL!
C
        IF((.NOT.NXT) .AND. GOPARR) THEN
           IPCOUNT = IPCOUNT + 1
           IF(MOD(IPCOUNT,NPROC).NE.ME) GO TO 580
        END IF
C
        J = KATOM(JJ)
        XJ = C(1,J)
        YJ = C(2,J)
        ZJ = C(3,J)
        J1 = KSTART(JJ)
        J2 = J1+KNG(JJ)-1
        LJT = KTYPE(JJ)
        MINJ = KMIN(JJ)
        MAXJ = KMAX(JJ)
        LOCJ = KLOC(JJ)-MINJ
        LJTDER=LJT+1
        NROOTS = (LITDER+LJTDER+NDER-2)/2 + 1
        RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
        IANDJ = II .EQ. JJ
C
C     ----- I PRIMITIVE
C
        DO 520 IG = I1,I2
        AI = EX(IG)
        ARRI = AI*RR
        AXI = AI*XI
        AYI = AI*YI
        AZI = AI*ZI
        CSI = CS(IG)
        CPI = CP(IG)
        CDI = CD(IG)
        CFI = CF(IG)
        CGI = CG(IG)
C
C     ----- J PRIMITIVE
C
        JGMAX = J2
        DO 500 JG = J1,JGMAX
          AJ = EX(JG)
          AA = AI+AJ
          AA1 = ONE/AA
          DUM = AJ*ARRI*AA1
          IF(DUM .GT. TOL) GO TO 500
          FAC = EXP(-DUM)
          CSJ = CS(JG)
          CPJ = CP(JG)
          CDJ = CD(JG)
          CFJ = CF(JG)
          CGJ = CG(JG)
          AX = (AXI+AJ*XJ)*AA1
          AY = (AYI+AJ*YJ)*AA1
          AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
          JMAX = MAXJ
          IJ = 0
          DUM1 = ZERO
          DUM2 = DUM1
          DO 360 I = MINI,MAXI
            IF(I.EQ.1)  DUM1=CSI*FAC
            IF(I.EQ.2)  DUM1=CPI*FAC
            IF(I.EQ.5)  DUM1=CDI*FAC
            IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.11) DUM1=CFI*FAC
            IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
            IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.21) DUM1=CGI*FAC
            IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
            IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
            IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
C
            IF(IANDJ) JMAX = I
            DO 360 J = MINJ,JMAX
              IF(J.EQ.1)  DUM2=DUM1*CSJ
              IF(J.EQ.2)  DUM2=DUM1*CPJ
              IF(J.EQ.5)  DUM2=DUM1*CDJ
              IF(J.EQ.8.AND.NORM) DUM2 = DUM2*SQRT3
              IF(J.EQ.11) DUM2=DUM1*CFJ
              IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
              IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
              IF(J.EQ.21) DUM2=DUM1*CGJ
              IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
              IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
              IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
C
              IJ = IJ+1
              IF(RMETHOD.EQ.RESC) THEN
                 DIJ(IJ)=DUM2
              ELSE
                 NDUM = IA(LOCI+I)+(LOCJ+J)
                 DEN=DENAB(NDUM)
                 IF(.NOT.IANDJ.OR.I.NE.J) DEN=DEN+DEN
                 DIJ(IJ)=FSC2/TWO*DUM2*DEN
C                THE ABOVE "TWO" COMES FROM NESC ITSELF
              ENDIF
C
C          THE DIFFERENCE IN FACTORS "2" BETWEEN THIS AND HELFEY
C          ROUTINES COMES FROM TRIANGULAR DENSITY ARRAY, OFF-DIAGONAL
C          PART OF WHICH IS DOUBLED IN HELFEY AND FOR RESC NOT DOUBLED
C          NOW, BUT LATER.
C
  360     CONTINUE
C
C     ..... HELLMANN-FEYNMAN TERM .....
C
          AAX = AA*AX
          AAY = AA*AY
          AAZ = AA*AZ
          DO 480 IC = 1,NAT
            ZNUC = -ZAN(IC)
            CX = C(1,IC)
            CY = C(2,IC)
            CZ = C(3,IC)
            XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
            IF(NROOTS.LE.3) CALL RT123
            IF(NROOTS.EQ.4) CALL ROOT4
            IF(NROOTS.EQ.5) CALL ROOT5
            IF(NROOTS.GE.6) CALL ROOT6
            DO 420 K = 1,NROOTS
              UU = AA*U(K)
              WW = W(K)*ZNUC
              WW=WW*(UU+UU)
              TT = ONE/(AA+UU)
              T = SQRT(TT)
              X0 = (AAX+UU*CX)*TT
              Y0 = (AAY+UU*CY)*TT
              Z0 = (AAZ+UU*CZ)*TT
              DO 400 J = 1,LJTDER
                NJ = J
                DO 400 I = 1,LITDER
                  NI = I
                  CALL VINT
                  XV(I,J,K) = XINT
                  YV(I,J,K) = YINT
                  ZV(I,J,K) = ZINT*WW
                  CALL DVINT
                  DXV(I,J,K) = XINT
                  DYV(I,J,K) = YINT
                  DZV(I,J,K) = ZINT*WW
  400         CONTINUE
              CALL QRDDERIJ( XVDIJ(1,1,K), YVDIJ(1,1,K), ZVDIJ(1,1,K),
     1                          XV(1,1,K),    YV(1,1,K),    ZV(1,1,K),
     2                      LIT,LJT,AI,AJ,5)
              CALL QRDDERIJ(DXVDIJ(1,1,K),DYVDIJ(1,1,K),DZVDIJ(1,1,K),
     1                         DXV(1,1,K),   DYV(1,1,K),   DZV(1,1,K),
     2                      LIT,LJT,AI,AJ,5)
  420       CONTINUE
           IJ=0
           DO 460 I=MINI,MAXI
            IX=IJX(I)
            IY=IJY(I)
            IZ=IJZ(I)
            JMAX=MAXJ
            IF(IANDJ) JMAX=I
            DO 460 J=MINJ,JMAX
             JX=IJX(J)
             JY=IJY(J)
             JZ=IJZ(J)
C
C     DUMIJI STANDS FOR PI*(DV/DRJ)*PI
C
             DUMXXX=ZERO
             DUMXYX=ZERO
             DUMXZX=ZERO
             DUMYXY=ZERO
             DUMYYY=ZERO
             DUMYZY=ZERO
             DUMZXZ=ZERO
             DUMZYZ=ZERO
             DUMZZZ=ZERO
             DO 430 K=1,NROOTS
C                            <PX (DV/DRX) PX> <V> <V>
               DUMXXX=DUMXXX+DXVDIJ(IX,JX,K)* YV(IY,JY,K)* ZV(IZ,JZ,K)
C                            <PX V PX> <DV/DRY> <V>
               DUMXYX=DUMXYX+ XVDIJ(IX,JX,K)*DYV(IY,JY,K)* ZV(IZ,JZ,K)
C                            <PX V PX> <V> <DV/DRZ> ET CETERA
               DUMXZX=DUMXZX+ XVDIJ(IX,JX,K)* YV(IY,JY,K)*DZV(IZ,JZ,K)
               DUMYXY=DUMYXY+DXV(IX,JX,K)* YVDIJ(IY,JY,K)* ZV(IZ,JZ,K)
               DUMYYY=DUMYYY+ XV(IX,JX,K)*DYVDIJ(IY,JY,K)* ZV(IZ,JZ,K)
               DUMYZY=DUMYZY+ XV(IX,JX,K)* YVDIJ(IY,JY,K)*DZV(IZ,JZ,K)
               DUMZXZ=DUMZXZ+DXV(IX,JX,K)* YV(IY,JY,K)* ZVDIJ(IZ,JZ,K)
               DUMZYZ=DUMZYZ+ XV(IX,JX,K)*DYV(IY,JY,K)* ZVDIJ(IZ,JZ,K)
               DUMZZZ=DUMZZZ+ XV(IX,JX,K)* YV(IY,JY,K)*DZVDIJ(IZ,JZ,K)
  430        CONTINUE
             IJ=IJ+1
C
C          I.E. <MU|P*(DV/DRJ)*P|NU>=SUM OVER I <MU|*PI*(DV/DRJ)*PI|NU>
C
             DUMX=DUMXXX+DUMYXY+DUMZXZ
             DUMY=DUMXYX+DUMYYY+DUMZYZ
             DUMZ=DUMXZX+DUMYZY+DUMZZZ
             DUM=DIJ(IJ)*AA1*PI212
             IF(RMETHOD.EQ.RESC) THEN
                NN=IA(LOCI+I)+(LOCJ+J)
                DVINTS(NN,1,IC)=DVINTS(NN,1,IC)+DUMX*DUM
                DVINTS(NN,2,IC)=DVINTS(NN,2,IC)+DUMY*DUM
                DVINTS(NN,3,IC)=DVINTS(NN,3,IC)+DUMZ*DUM
             ELSE
               DE(1,IC) = DE(1,IC)+DUM*DUMX
               DE(2,IC) = DE(2,IC)+DUM*DUMY
               DE(3,IC) = DE(3,IC)+DUM*DUMZ
             ENDIF
  460       CONTINUE
  480     CONTINUE
  500 CONTINUE
  520 CONTINUE
C
C     ----- END OF *PRIMITIVE* LOOPS -----
C
  580 CONTINUE
  600 CONTINUE
C
C     ----- END OF *SHELL* LOOPS -----
C
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
C
C     THE RESULTS ARE NOT SAVED UNTIL THE AO DERIVATIVE CONTRIBUTIONS
C     ARE ADDED
C
      RETURN
C
 9000 FORMAT(/10X,22(1H-)/10X,'HELLMANN-FEYNMAN PVP FORCE'/10X,22(1H-))
      END
C*MODULE QREL     *DECK GDERI
      SUBROUTINE GDERI(DXDI,DYDI,DZDI,X,Y,Z,LD,LIT,LJT,AI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION  X  (LD,1), Y  (LD,1), Z  (LD,1)
      DIMENSION DXDI(5,1),DYDI(5,1),DZDI(5,1)
C
C     A CLONE OF DERI WITH LD ADDED
C
C     ----- COMPUTE DERIVATIVE INTEGRALS <MU'|V|NU>
C
      DO 110 J=1,LJT
        DXDI(1,J)=X(2,J)*(AI+AI)
        DYDI(1,J)=Y(2,J)*(AI+AI)
        DZDI(1,J)=Z(2,J)*(AI+AI)
  110 CONTINUE
C
      IF(LIT.EQ.1) RETURN
C
      DO 120 I=2,LIT
        DO 120 J=1,LJT
          DXDI(I,J)=X(I+1,J)*(AI+AI) - X(I-1,J)*(I-1)
          DYDI(I,J)=Y(I+1,J)*(AI+AI) - Y(I-1,J)*(I-1)
          DZDI(I,J)=Z(I+1,J)*(AI+AI) - Z(I-1,J)*(I-1)
  120 CONTINUE
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRELX
      SUBROUTINE QRELX(HAO,SAO,TAO,PVP,L1,L2,DBUG,ISAVE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      DIMENSION HAO(*),SAO(*),TAO(*),PVP(*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /FMCOM / X(1)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DATA TWO /2.0D+00/
      DATA RESC /8HRESC    /
C
C     ----- DOUGLAS-KROLL APPROACH -----
C     ----- RELATIVISTIC SCHEME BY ELIMINATING SMALL-COMPONENTS -----
C
C           CODED BY TAKAHITO NAKAJIMA
C
C           DEPARTMENT OF APPLIED CHEMISTRY, UNIVERSITY OF TOKYO (2001).
C     CHECK THE GRADIENT CALCULATION
      CALL DERCHK(NDER)
      L3 = L1 * L1
      CALL VALFM(LOADFM)
      LHQR    = LOADFM  + 1
      LVAO    = LHQR    + L3
      LTAO2   = LVAO    + L2
      LXX     = LTAO2   + L2
      LSDIAG  = LXX     + L3
      LUU     = LSDIAG  + L1
      LTDIAG  = LUU     + L3
      LEP     = LTDIAG  + L1
      LAP     = LEP     + L1
      LRP     = LAP     + L1
      LXU     = LRP     + L1
      LSXU    = LXU     + L3
      LEVEN1  = LSXU    + L3
      LRE1R   = LEVEN1  + L3
      LW1SQ   = LRE1R   + L3
      LTWRK0  = LW1SQ   + L3
      LTWRK1A = LTWRK0  + L1
      LTWRK2A = LTWRK1A + L2
      LTWRK2B = LTWRK2A + L3
      LTWRK2C = LTWRK2B + L3
      LGIVEB  = LTWRK2C + L3
      LINDB   = LGIVEB  + L1 * 8
      LAST    = LINDB + (L1 - 1) / NWDVAR + 1
      NEED    = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- SPECIFY RELATIVISTIC METHOD -----
C
      IMETHOD=IQRORD
      IF (RMETHOD .EQ. RESC)    IMETHOD = 0
      KAMERAD=MOD(MODQR/2,2)
C
C     ----- CONSTRUCT VAO AND TAO*2 -----
C
      CALL VSUB(TAO,1,HAO,1,X(LVAO),1,L2)
      CALL DCOPY(L2,TAO,1,X(LTAO2),1)
      CALL DSCAL(L2,TWO,X(LTAO2),1)
      IF (NDER .GT. 0) THEN   ! FOR GRADIENT CALCULATION
         CALL DAWRIT(IDAF,IODA,X(LTAO2),L2,NDARELB+4,0)
C        CALL DAWRIT(IDAF,IODA,PVP,L2,NDARELB+5,0)
         CALL DAWRIT(IDAF,IODA,X(LVAO),L2,NDARELB+6,0)
C        CALL DAWRIT(IDAF,IODA,SAO,L2,NDARELB+20,0)
      END IF
      IF (DBUG) THEN
         WRITE(IW,*)'$$$$$ DARWIN INTEGRALS $$$$$',L1,L2
         CALL PRTRI(PVP,L1)
      ENDIF
C
C     ----- CALCULATE RELATIVISTIC HAMILTONIAN -----
C
      CALL HQREL(HAO,SAO,X(LTAO2),X(LVAO),PVP,X(LHQR),X(LXX),X(LSDIAG),
     *           X(LUU),X(LTDIAG),X(LEP),X(LAP),X(LRP),X(LXU),X(LSXU),
     *           X(LEVEN1),X(LRE1R),X(LW1SQ),
     *           X(LTWRK0),X(LTWRK1A),X(LTWRK2A),X(LTWRK2B),X(LTWRK2C),
     *           X(LGIVEB),X(LINDB),NDER,L1,L2,KAMERAD,IMETHOD)
C
C     ----- WRITE RELATIVISTIC HAMILTONIAN TO DAF -----
C
      CALL DAWRIT(IDAF,IODA,HAO,L2,11,0)
      IF (ISAVE .EQ. 1) CALL DAWRIT(IDAF,IODA,HAO,L2,311,0)
      IF (DBUG) THEN
         WRITE(IW,*)'+++++ RELATIVISTIC HAO +++++',CLIG,CLIG2
         CALL PRTRIL(HAO,L1)
      END IF
C
      CALL RETFM(NEED)
C
      RETURN
      END
C
C*MODULE QREL     *DECK HQREL
      SUBROUTINE HQREL(HAO,SAO,TAO2,VAO,PVP,HQR,XX,SDIAG,UU,TDIAG,
     *                 EP,AP,RP,XU,SXU,EVEN1,RE1R,W1SQ,
     *                 TWRK0,TWRK1A,TWRK2A,TWRK2B,TWRK2C,
     *                 GIVEB,INDB,NDER,L1,L2,KAMERAD,IMETHOD)
C
C     ----- DOUGLAS-KROLL APPROACH -----
C     ----- RELATIVISTIC SCHEME BY ELIMINATING SMALL-COMPONENTS -----
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION HAO(*),SAO(*),TAO2(*),VAO(*),PVP(*),HQR(L1,*)
      DIMENSION XX(L1,*),SDIAG(*),UU(L1,*),TDIAG(*)
      DIMENSION AP(*),RP(*),EP(*),XU(L1,*),SXU(L1,*)
      DIMENSION EVEN1(L1,*),RE1R(L1,*),W1SQ(L1,*)
      DIMENSION TWRK0(*),TWRK1A(*)
      DIMENSION TWRK2A(L1,*),TWRK2B(L1,*),TWRK2C(L1,*)
      DIMENSION GIVEB(L1,8),INDB(L1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DATA HALF/0.5D+00/,ONE/1.0D+00/
C
      CALL QRSXU(SXU,SAO,TAO2,XX,SDIAG,UU,TDIAG,XU,TWRK0,TWRK1A,GIVEB,
     $           INDB,NDER,KAMERAD,L1,L0,L2)
      L02=(L0*L0+L0)/2
C
C     ----- CONSTRUCT RELATIVISTIC HAMILTONIAN -----
C
      CALL VCLR(HQR,1,L1*L1)
C
C     ----- CALCULATE SQUARE-ROOT (SQR) TERM -----
C
      DO I = 1, L0
         ID=(I-1)*(L1-L0)
         EP(I) = CLIG * SQRT(TDIAG(I) + CLIG2)
         HQR(I-ID,I) = EP(I) - CLIG2
C        HQR(I,I) = EP(I) - CLIG2
C        THE UGLY USE OF -(I-1)*(L1-L0) IS TO MAKE HQR TO BE (L0,L0)
C        DESPITE DECLARING IT AS (L1,L1)
      END DO
C
C     ----- CALCULATE KINEMATICAL FACTORS -----
C
C     DK
      IF (IMETHOD .GE. 1) THEN
         DO I = 1, L0
            AP(I) = SQRT((EP(I) + CLIG2) * HALF / EP(I))
            RP(I) = CLIG / (EP(I) + CLIG2)
         END DO
      ELSE IF (IMETHOD .EQ. 0) THEN
C     RESC: SECRET OPTION $RELWFN NORDER=0
         DO I = 1, L0
            AP(I) = ( SQRT(ONE + CLIG2 * TDIAG(I)
     *            / ((EP(I) + CLIG2) * (EP(I) + CLIG2))))
     *            / (EP(I) + CLIG2)
C           O OPERATOR
            RP(I) = CLIG / (EP(I) + CLIG2)
C           Q OPERATOR
         END DO
      ELSE
         CALL ABRT
      END IF
C
C     ----- CONSTRACT (SXU) APRP (XU)T FOR SPIN-ORBIT CALCULATION -----
C
      CALL VCLR(TWRK2A,1,L1*L0)
      DO I = 1, L0
         TWRK2A(I,I) = AP(I) * RP(I)
      END DO
      CALL MRARBR(SXU,L1,L1,L0,TWRK2A,L1,L0,TWRK2B,L1)
      CALL MRARTR(TWRK2B,L1,L1,L0,XU,L1,L1,TWRK2A,L1)
      NDARE=NDARELB
      CALL DAWRIT(IDAF,IODA,TWRK2A,L1*L1,NDARE  ,0)
      CALL DAWRIT(IDAF,IODA,TWRK2A,L1*L1,NDARE+1,0)
C
C   ----- SAVE RELATIVISTIC KINETIC ENERGY (SXU)EI*(SXU)T
C
C     CALL DAREAD(IDAF,IODA,TWRK1A,L2,13,0)
C     WRITE(6,*) 'READING REL KIN',TWRK1A(1),TWRK1A(2)
      CALL MRARBR(SXU,L1,L1,L0,HQR,L0,L0,TWRK2A,L1)
      CALL MRARTR(TWRK2A,L1,L1,L0,SXU,L1,L1,TWRK2B,L1)
C     WRITE(6,*) 'HAVING HQR',HQR(1,1),TWRK2B(1,1),SXU(1,1)
      CALL CPYSQT(TWRK2B,TWRK1A,L1,1)
C     WRITE(6,*) 'SAVING REL KIN',TWRK1A(1),TWRK1A(2)
      CALL DAWRIT(IDAF,IODA,TWRK1A,L2,13,0)
C
C     ----- TRANSFORM VAO AND PVP TO P-SPACE BASIS -----
C
      CALL TFTRI(TWRK2A,VAO,XU,TWRK0,L0,L1,L1)
      CALL DCOPY(L02,TWRK2A,1,VAO,1)
C
      CALL TFTRI(TWRK2A,PVP,XU,TWRK0,L0,L1,L1)
      CALL DCOPY(L02,TWRK2A,1,PVP,1)
C
C     ----- CALCULATE DK1 OR RESC -----
C
      CALL DK1X(HQR,VAO,PVP,AP,RP,EVEN1,RE1R,TDIAG,L0,CLIG,IMETHOD)
      IF (IMETHOD .EQ. 0 .OR. IMETHOD .EQ. 1) GO TO 1000
C
C     ----- DEVIDE VAO AND PVP (IN P-SPACE BASIS) BY EP + EP' -----
C           FOR LATER CALCULATIONS.
C
      IJ = 0
      DO I = 1, L0
         DO J = 1, I
            IJ = IJ + 1
            TEMP = ONE / (EP(I) + EP(J))
            VAO(IJ) = VAO(IJ) * TEMP
            PVP(IJ) = PVP(IJ) * TEMP
         END DO
      END DO
C
C     ----- CALCULATE DK2 AND DK2FULL -----
C
      CALL DK2X(HQR,VAO,PVP,EP,AP,RP,TDIAG,W1SQ,
     *          TWRK2A,TWRK2B,L0,IMETHOD)
      IF (IMETHOD .EQ. 2 .OR. IMETHOD .EQ. 12) GO TO 1000
C
C     ----- CALCULATE DK3 -----
C
      CALL DK3X(HQR,VAO,PVP,AP,RP,TDIAG,EVEN1,RE1R,W1SQ,
     *          TWRK2A,TWRK2B,TWRK2C,L0,IMETHOD)
      IF (IMETHOD .EQ. 3 .OR. IMETHOD .EQ. 13) GO TO 1000
C
C     ----- BACK-TRANSFORM TO R-SPACE BASIS -----
C           BY MULTIPLING SXU
C
 1000 CONTINUE
      CALL MRARBR(SXU,L1,L1,L0,HQR,L0,L0,TWRK2A,L1)
      CALL MRARTR(TWRK2A,L1,L1,L0,SXU,L1,L1,TWRK2B,L1)
      CALL SYMTRZE(TWRK2B,L1,L1)
      CALL CPYSQT(TWRK2B,HAO,L1,1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK DK1X
      SUBROUTINE DK1X(HQR,VAO,PVP,AP,RP,EVEN1,RE1R,TDIAG,
     *                L1,CLIG,IMETHOD)
C
C     ----- DK1 HAMILTONIAN (FREE-PARTICLE) -----
C     ----- RELATIVISTIC SCHEME BY ELIMINATING SMALL-COMPONENTS -----
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAO(*),PVP(*),HQR(L1,*)
      DIMENSION AP(*),RP(*),TDIAG(*)
      DIMENSION EVEN1(L1,*),RE1R(L1,*)
      DATA HALF/0.5D+00/,ONE/1.0D+00/
C
C     ----- CALCULATE A V A -----
C           ALSO, CALCULATE AR^2 V R^2A = A K^2 P^2 V P^2 K^2 A
C           FOR HIGHER-ORDER DK
C
      IF (IMETHOD .GE. 1) THEN
C     DK
         IJ = 0
         DO I = 1, L1
            DO J = 1, I
               IJ = IJ + 1
               EVEN1(I,J) = AP(I) * VAO(IJ) * AP(J)
               RE1R(I,J)  = RP(I) * RP(I) * TDIAG(I) * EVEN1(I,J)
     *                    * RP(J) * RP(J) * TDIAG(J)
            END DO
         END DO
      ELSE IF (IMETHOD .EQ. 0) THEN
C     RESC
         IJ = 0
         DO I = 1, L1
            DO J = 1, I
               IJ = IJ + 1
               EVEN1(I,J) = CLIG * VAO(IJ) *
     *                    ( SQRT(RP(I) * RP(J)) * AP(I) / AP(J)
     *                    + SQRT(RP(J) * RP(I)) * AP(J) / AP(I))
            END DO
         END DO
      END IF
C
C     ----- CALCULATE AR V RA -----
C
      IF (IMETHOD .GE. 1) THEN
C     DK
         IJ = 0
         DO I = 1, L1
            DO J = 1, I
               IJ = IJ + 1
               EVEN1(I,J) = EVEN1(I,J)
     *                    + AP(I) * RP(I) * PVP(IJ) * RP(J) * AP(J)
               EVEN1(J,I) = EVEN1(I,J)
               RE1R(I,J)  = RE1R(I,J)
     *                    + AP(I) * RP(I) * PVP(IJ) * RP(J) * AP(J)
               RE1R(J,I)  = RE1R(I,J)
            END DO
         END DO
      ELSE IF (IMETHOD .EQ. 0) THEN
C     RESC
         IJ = 0
         DO I = 1, L1
            DO J = 1, I
               IJ = IJ + 1
               EVEN1(I,J) = EVEN1(I,J)
     *                    + PVP(IJ) * (RP(I) * RP(J) * AP(I) / AP(J)
     *                               + RP(J) * RP(I) * AP(J) / AP(I))
     *                    * HALF
               EVEN1(J,I) = EVEN1(I,J)
            END DO
         END DO
      END IF
C
C     ----- ADD EVEN1 TO HQR -----
C
      CALL DAXPY(L1*L1,ONE,EVEN1,1,HQR,1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK DK2X
      SUBROUTINE DK2X(HQR,VAO,PVP,EP,AP,RP,TDIAG,W1SQ,
     *                TWRK2A,TWRK2B,L1,IMETHOD)
C
C     ----- DK2 HAMILTONIAN -----
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAO(*),PVP(*),HQR(L1,*)
      DIMENSION EP(*),AP(*),RP(*),TDIAG(*)
      DIMENSION W1SQ(L1,*)
      DIMENSION TWRK2A(L1,*),TWRK2B(L1,*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DATA HALF/0.5D+00/
C
C     ----- CALCULATE - 1/2 (EPW1W1 + W1W1EP) -----
C
C     - EPW1W1 = - EP AK(PV'P)KA AV'A
C                - EP AV'A AK(PV'P)KA
C                + EP AV'A P^2K^2 AV'A
C                + EP AK(PV'P)KA P^-2K^-2 AK(PV'P)KA
C                - EP AK(PXV'P)KA P^-2K^-2 AK(PXV'P)KA
C
C     - W1W1EP = - AK(PV'P)KA AV'A EP
C                - AV'A AK(PV'P)KA EP
C                + AV'A P^2K^2 AV'A EP
C                + AK(PV'P)KA P^-2K^-2 AK(PV'P)KA EP
C                - AK(PXV'P)KA P^-2K^-2 AK(PXV'P)KA EP
C
C     ------- CALCULATE W1W1
C
C     ------- TERM 1 IN W1W1
C
C     AK(PV'P)KA AV'A
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * RP(I) * PVP(IJ) * RP(J) * AP(J)
            TWRK2A(J,I) = TWRK2A(I,J)
            TWRK2B(I,J) = AP(I) * VAO(IJ) * AP(J)
            TWRK2B(J,I) = TWRK2B(I,J)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2A,L1,TWRK2B,L1,0.0D+00,
     *           W1SQ,L1)
C
C     ------- TERM 2 IN W1W1
C
C     AV'A AK(PV'P)KA
C
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2B,L1,TWRK2A,L1,1.0D+00,
     *           W1SQ,L1)
C
C     ------- TERM 3 IN W1W1
C
C     - AV'A P^2K^2 AV'A
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * VAO(IJ)
     *                  * AP(J) * TDIAG(J) * RP(J) * RP(J)
            TWRK2A(J,I) = AP(J) * VAO(IJ)
     *                  * AP(I) * TDIAG(I) * RP(I) * RP(I)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,-1.0D+00,TWRK2A,L1,TWRK2B,L1,1.0D+00,
     *           W1SQ,L1)
C
C     ------- TERM 4 IN W1W1
C
C     - AK(PV'P)KA P^-2K^-2 AK(PV'P)KA
C     = - AK(PV'P)A P^-2 A(PV'P)KA
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * RP(I) * PVP(IJ) * AP(J) / TDIAG(J)
            TWRK2A(J,I) = AP(J) * RP(J) * PVP(IJ) * AP(I) / TDIAG(I)
            TWRK2B(I,J) = AP(I) * PVP(IJ) * RP(J) * AP(J)
            TWRK2B(J,I) = AP(J) * PVP(IJ) * RP(I) * AP(I)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,-1.0D+00,TWRK2A,L1,TWRK2B,L1,1.0D+00,
     *           W1SQ,L1)
C
C     ------- TERM 5 IN W1W1
C             CROSS-PRODUCT TERM
C
      IF (IMETHOD .GE. 12) THEN
         WRITE(IW,*)'+++++ CURRENTLY, NOT YET IMPLEMENTED +++++'
      END IF
CNAKA      WRITE(IW,*)'+++++ W1SQ IN DK2X +++++'
CNAKA      CALL MATOUT(W1SQ,L1,L1)
C
C     ------- CALCULATE - 1/2 (EP * W1W1 + W1W1 * EP)
C
      DO I = 1, L1
         DO J = 1, I
            HQR(I,J) = HQR(I,J)
     *               - (EP(I) * W1SQ(I,J) + W1SQ(I,J) * EP(J)) * HALF
            HQR(J,I) = HQR(I,J)
         END DO
      END DO
CNAKA      WRITE(IW,*)'+++++ HQR IN DK2X +++++'
CNAKA      CALL MATOUT(HQR,L1,L1)
C
C     ----- CALCULATE - W1EPW1 -----
C
C     - W1EPW1 = - AK(PV'P)KA EP AV'A
C                - AV'A EP AK(PV'P)KA
C                + AV'A P^2K^2 EP AV'A
C                + AK(PV'P)KA EP K^-2P^-2 AK(PV'P)KA
C                - AK(PXV'P)KA EP K^-2P^-2 AK(PXV'P)KA
C
C     ------- TERM 1 IN - W1EPW1
C
C     - AK(PV'P)KA EP AV'A
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * RP(I) * PVP(IJ)
     *                  * RP(J) * AP(J) * EP(J)
            TWRK2A(J,I) = AP(J) * RP(J) * PVP(IJ)
     *                  * RP(I) * AP(I) * EP(I)
            TWRK2B(I,J) = AP(I) * VAO(IJ) * AP(J)
            TWRK2B(J,I) = TWRK2B(I,J)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,-1.0D+00,TWRK2A,L1,TWRK2B,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 2 IN - W1EPW1
C
C     - AV'A EP AK(PV'P)KA
C
      CALL DGEMM('N','T',L1,L1,L1,-1.0D+00,TWRK2B,L1,TWRK2A,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 3 IN - W1EPW1
C
C     + AV'A P^2K^2 EP AV'A
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * VAO(IJ)
     *                  * AP(J) * TDIAG(J) * RP(J) * RP(J) * EP(J)
            TWRK2A(J,I) = AP(J) * VAO(IJ)
     *                  * AP(I) * TDIAG(I) * RP(I) * RP(I) * EP(I)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2A,L1,TWRK2B,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 4 IN - W1EPW1
C
C     + AK(PV'P)KA EP K^-2P^-2 AK(PV'P)KA
C     = + AK(PV'P)A EP P^-2 A(PV'P)KA
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * RP(I) * PVP(IJ)
     *                  * AP(J) * EP(J) / TDIAG(J)
            TWRK2A(J,I) = AP(J) * RP(J) * PVP(IJ)
     *                  * AP(I) * EP(I) / TDIAG(I)
            TWRK2B(I,J) = AP(I) * PVP(IJ) * RP(J) * AP(J)
            TWRK2B(J,I) = AP(J) * PVP(IJ) * RP(I) * AP(I)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2A,L1,TWRK2B,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 5 IN - W1EPW1
C             CROSS-PRODUCT TERM
C
      IF (IMETHOD .GE. 12) THEN
         WRITE(IW,*)'+++++ CURRENTLY, NOT YET IMPLEMENTED +++++'
      END IF
CNAKA      WRITE(IW,*)'+++++ HQR IN DK2X +++++'
CNAKA      CALL MATOUT(HQR,L1,L1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK DK3X
      SUBROUTINE DK3X(HQR,VAO,PVP,AP,RP,TDIAG,EVEN1,RE1R,W1SQ,
     *                TWRK2A,TWRK2B,TWRK2C,L1,IMETHOD)
C
C     ----- DK3 HAMILTONIAN -----
C
C     HDK3 = HDK2 + 1/2 [W1,[W1,E1]]
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAO(*),PVP(*),HQR(L1,*)
      DIMENSION AP(*),RP(*),TDIAG(*)
      DIMENSION EVEN1(L1,*),RE1R(L1,*)
      DIMENSION TWRK2A(L1,*),TWRK2B(L1,*),TWRK2C(L1,*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- CALCULATE + 1/2 (E1W1W1 + W1W1E1) -----
C
C     + E1W1W1 = + E1 AK(PV'P)KA AV'A
C                + E1 AV'A AK(PV'P)KA
C                - E1 AV'A P^2K^2 AV'A
C                - E1 AK(PV'P)KA K^-2P^-2 AK(PV'P)KA
C                - AK(PXVP)KA AK(PXV'P)KA AV'A                  ! TERM 1
C                - AK(PXVP)KA AV'A AK(PXV'P)KA                  ! TERM 2
C                + E1 AK(PXV'P)KA K^-2P^-2 AK(PXV'P)KA          ! TERM 4
C                + AK(PXVP)KA AK(PXV'P)KA K^-2P^-2 AK(PV'P)KA   ! TERM 4
C                + AK(PXVP)KA AK(PV'P)KA K^-2P^-2 AK(PXV'P)KA   ! TERM 4
C              = + E1 W1^2
C                - AK(PXVP)KA AK(PXV'P)KA AV'A                  ! TERM 1
C                - AK(PXVP)KA AV'A AK(PXV'P)KA                  ! TERM 2
C                + AK(PXVP)KA AK(PXV'P)KA K^-2P^-2 AK(PV'P)KA   ! TERM 4
C                + AK(PXVP)KA AK(PV'P)KA K^-2P^-2 AK(PXV'P)KA   ! TERM 4
C
C     + W1W1E1 = + AK(PV'P)KA AV'A E1
C                + AV'A AK(PV'P)KA E1
C                - AV'A P^2K^2 AV'A E1
C                - AK(PV'P)KA K^-2P^-2 AK(PV'P)KA E1
C                - AV'A AK(PXV'P)KA AK(PXVP)KA                  ! TERM 2
C                - AK(PXV'P)KA AV'A AK(PXVP)KA                  ! TERM 1
C                + AK(PXV'P)KA K^-2P^-2 AK(PXV'P)KA E1          ! TERM 4
C                + AK(PV'P)KA K^-2P^-2 AK(PXV'P)KA AK(PXVP)KA   ! TERM 4
C                + AK(PXV'P)KA K^-2P^-2 AK(PV'P)KA AK(PXVP)KA   ! TERM 4
C              = + W1^2 E1
C                - AV'A AK(PXV'P)KA AK(PXVP)KA                  ! TERM 2
C                - AK(PXV'P)KA AV'A AK(PXVP)KA                  ! TERM 1
C                + AK(PV'P)KA K^-2P^-2 AK(PXV'P)KA AK(PXVP)KA   ! TERM 4
C                + AK(PXV'P)KA K^-2P^-2 AK(PV'P)KA AK(PXVP)KA   ! TERM 4
C
C     ------- CALCULATE + 1/2 (E1 * W1W1 + W1W1 * E1)
C
      CALL DGEMM('N','N',L1,L1,L1,0.5D+00,EVEN1,L1,W1SQ,L1,
     *           1.0D+00,HQR,L1)
      CALL DGEMM('N','N',L1,L1,L1,0.5D+00,W1SQ,L1,EVEN1,L1,
     *           1.0D+00,HQR,L1)
C
C     ------- CROSS-PRODUCT TERMS IN + 1/2 (E1W1W1 + W1W1E1)
C
      IF (IMETHOD .GE. 13) THEN
         WRITE(IW,*)'+++++ CURRENTLY, NOT YET IMPLEMENTED +++++'
      END IF
C
C     ----- CALCULATE - W1E1W1 -----
C
C     - W1E1W1 = - AK(PV'P)KA R^-2 RE(1)R AV'A
C                - AV'A RE1R R^-2 AK(PV'P)KA
C                + AV'A RE1R AV'A
C                + AK(PV'P)KA R^-2 RE1R R^-2 AK(PV'P)KA
C                + AK(PXV'P)KA R^-2 AK(PXVP)KA AVG_AA           ! TERM 1
C                + AVG_AA AK(PXVP)KA R^-2 K(PXV'P)KA            ! TERM 2
C                - AK(PXV'P)KA R^-2 RE1R R^-2 AK(PXV'P)KA       ! TERM 4
C                - AK(PXV'P)KA R^-2 AK(PXVP)KA R^-2 AKVDG_AKA   ! TERM 4
C                - AKVDG_AKA R^-2 AK(PXVP)KA R^-2 AK(PXV'P)KA   ! TERM 4
C
C     ------- TERM 1 IN - W1E1W1
C
C     - AK(PV'P)KA R^-2 RE(1)R AV'A
C     = - KA(PV'P)AK^-1P^-2 RE1R AV'A
C
      IJ = 0
      DO I = 1, L1
         DO J = 1, I
            IJ = IJ + 1
            TWRK2A(I,J) = AP(I) * RP(I) * PVP(IJ)
     *                  * AP(J) / (RP(J) * TDIAG(J))
            TWRK2A(J,I) = AP(J) * RP(J) * PVP(IJ)
     *                  * AP(I) / (RP(I) * TDIAG(I))
            TWRK2B(I,J) = AP(I) * VAO(IJ) * AP(J)
            TWRK2B(J,I) = TWRK2B(I,J)
         END DO
      END DO
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2A,L1,RE1R,L1,0.0D+00,
     *           TWRK2C,L1)
      CALL DGEMM('N','N',L1,L1,L1,-1.0D+00,TWRK2C,L1,TWRK2B,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 2 IN - W1E1W1
C
C     - AV'A RE1R R^-2 AK(PV'P)KA
C     = - AV'A RE1R P^-2K^-1A(PV'P)AK
C
      CALL DGEMM('T','T',L1,L1,L1,-1.0D+00,TWRK2B,L1,TWRK2C,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 3 IN - W1E1W1
C
C     + AV'A RE1R AV'A
C
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2B,L1,RE1R,L1,0.0D+00,
     *           TWRK2C,L1)
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2C,L1,TWRK2B,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- TERM 4 IN - W1E1W1
C
C     + AK(PV'P)KA R^-2 RE1R R^-2 AK(PV'P)KA
C     = + KA(PV'P)A K^-1P^-2 RE1R P^-2K^-1A(PV'P)AK
C
      CALL DGEMM('N','N',L1,L1,L1,1.0D+00,TWRK2A,L1,RE1R,L1,0.0D+00,
     *           TWRK2C,L1)
      CALL DGEMM('N','T',L1,L1,L1,1.0D+00,TWRK2C,L1,TWRK2A,L1,1.0D+00,
     *           HQR,L1)
C
C     ------- CROSS-PRODUCT TERMS IN - W1E1W1
C
      IF (IMETHOD .GE. 13) THEN
         WRITE(IW,*)'+++++ CURRENTLY, NOT YET IMPLEMENTED +++++'
      END IF
CNAKA      WRITE(IW,*)'+++++ HQR IN DK3X +++++'
CNAKA      CALL MATOUT(HQR,L1,L1)
C
      RETURN
      END
C
C*MODULE QREL     *DECK QRNUMDER
      SUBROUTINE QRNUMDER(DENAB,L2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXSH=5000, MXGTOT=20000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (NIODA=950)
      PARAMETER (HALF=0.5D+00,ONE=1.0D+00)
C
      LOGICAL UNCON,DBG,GOPARR,DSKWRK,MASWRK,
     *        SVDSKW,SVGPAR,NXT,SKIP,CART,SOME
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /CONV  / ACURCY,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /DAIOLN/ IRECLN,IRECST,IFILEN(NIODA)
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(NIODA)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /MASSES/ ZMASS(MXATM)
C
      DIMENSION DENAB(*),SSIGN(2)
      DATA SSIGN/1.0D+00,-1.0D+00/
      DATA DEBUG/8HDEBUG   /, DBUGME/8HQRNUMDER/,DBG1/8HGRD1    /
C
C     ----- MEMORY ALLOCATION -----
C
      DBG = EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME .OR. EXETYP.EQ.DBG1
C     THIS TOL IS ONLY USED FOR LINEAR DEPENDENCY
      TOL=1.0D-04
      UNCON = MOD(MODQR,2).EQ.1
      KAMERAD=MOD(MODQR/2,2)
      CART=MOD(MODQR/131072,2).EQ.1
C
      IF(KAMERAD.NE.0.AND.ISPHER.EQ.1.AND.MASWRK) WRITE(IW,8000)
C
C        CALL ABRT
C     THE REASON FOR THE HONDO OPTION FAILURE IS NOT TOTALLY CLEAR.
C     THE ENERGY MATCHES NWCHEM, SO THE PROBLEM IS PROBABLY LIMITED
C     TO THE GRADIENT. PERHAPS NWCHEM DOES NOT USE Q-MATRIX AND THE
C     Q-MATRIX EFFECT (Q-MATRIX DEPENDENCE ON THE GEOMETRY IS THE
C     CAUSE OF EVIL?
C     THE CODE WILL RUN, HOWEVER, AND PRODUCE SOME GRADIENT.
C
C     CONSTRUCT TOTALLY SYMMETRIC OFFSETS
C     AT PRESENT FRAGMENTS ARE NOT SUPPORTED
C
      IF(NFRG.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'QRNUMDER: NO EFP CODE EXISTS'
         CALL ABRT
      END IF
C
      NCCF = 3*NAT + 6*NFRG
      CALL VALFM(LOADFM)
      LX   = LOADFM+ 1
      LGTOT= LX    + NCCF
      LRM  = LGTOT + NCCF
      LP   = LRM   + 3*NAT + 21*NFRG
      LGE  = LP    + NCCF*NCCF
      LDR  = LGE   + 3*MXATM
      LT   = LDR   + NCCF
      LWRK = LT    + NCCF*NCCF
      LAST = LWRK  + NCCF
      NEEDP= LAST-LOADFM-1
      CALL GETFM(NEEDP)
      IF(CART) THEN
C       CARTESIAN SHIFT MODE. P IS A UNIT MATRIX.
        CALL VCLR(X(LP),1,NCCF*NCCF)
        DO I=1,NCCF
           X(LP+I-1+(I-1)*NCCF)=ONE
        ENDDO
        NMODE=NCCF
      ELSE
        CALL VCLR(X(LGE),1,3*MXATM)
C       FROM PRJGRD WE ONLY WANT THE P MATRIX (AND THE REDUCED MASSES).
        CALL PRJGRD(X(LGE),C,X(LX),X(LGTOT),X(LRM),X(LP),ZMASS,NAT,
     *              3*NAT,NCCF)
C       UN-MASS-WEIGH
        DO I=1,NCCF
          CALL DSCAL(NCCF,ONE/X(LRM+I-1),X(LP+(I-1)*NCCF),1)
        ENDDO
C       PROJECT EACH ROW ON THE TOTALLY SYMMETRIC REPRESENTATION
        DO I=1,NCCF
          CALL DCOPY(NCCF,X(LP+I-1),NCCF,X(LDR),1)
          CALL SYMDR(X(LDR))
          RABS=SQRT(DDOT(NCCF,X(LDR),1,X(LDR),1))
C         NORMALISE
          IF(RABS.GT.TOL) CALL DSCAL(NCCF,ONE/RABS,X(LDR),1)
          CALL DCOPY(NCCF,X(LDR),1,X(LP+I-1),NCCF)
        ENDDO
C
C       THROW AWAY TRANSLATIONS, ROTATIONS AND NONSYMMETRIC STRETCHES.
C       STRETCHES ARE STORED AS ROWS, SO TRANSPOSE NOW AND BACK AFTER.
C       (ROWS BECAUSE COORDINATES ARE TRANSFORMED X'=AX: ROWS OF A;
C       GRADIENT IS D/DX=D/DX' * DX'/DX = D/DX' * A).
C
        CALL TRPOSQ(X(LP),NCCF)
        CALL UNLINDEP(NCCF,NCCF,X(LP),NCCF,X(LT),NCCF,X(LWRK),TOL,NMODE)
        CALL TRPOSQ(X(LP),NCCF)
C       WRITE(6,*) 'FINAL P MATRIX IS',NCCF,NMODE
C       CALL PRSQ(X(LP),NCCF,NMODE,NCCF)
      ENDIF
C
      L1 = NUM
      LU2 = (NUMU * NUMU + NUMU) / 2
      L3 = L1*L1
      NVIB=2
      IF(MOD(MODQR/65536,2).EQ.1) NVIB=1
      NQRI=3-NVIB
      NBUF1=0
      NBUF2=0
      NBUF3=0
      CALL VALFM(LOADFM)
      LHR=LOADFM+1
      LHNR=   LHR+L2
      LHNRUNC=LHNR+L2
      LTAO2=  LHNRUNC+LU2
      LQRINT= LTAO2+LU2
      LDQRINT=LQRINT+L2*NQRI
      LDE=    LDQRINT+L2
      LDEC=   LDE+3*NAT
      LAST=   LDEC+3*NAT
C     IF(GOPARR) THEN
        IF(UNCON) THEN
          NBUF1=MXGTOT*6
          NBUF2=(MXSH*7+1)/NWDVAR+1
          NBUF3=NUMU*NUM
C         NBUF3 IS ALSO USED FOR NUM*NUM
        ENDIF
        IF(KAMERAD.EQ.1.AND..NOT.UNCON) NBUF3=L3
        LIUDA=  LAST
C       NIODA=950 - SAME NUMBER AS IN IOFILE, IODA.
        LJUDAS= LIUDA+(NIODA-1)/NWDVAR+1
        LFILEN= LJUDAS+(NIODA-1)/NWDVAR+1
        LBUF=   LFILEN+(NIODA-1)/NWDVAR+1
        LAST=   LBUF+MAX(NBUF1,NBUF2,NBUF3)
C     ENDIF
      NEED=   LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      SOME = NPRINT.NE.-5  .AND.  MASWRK
      NPRINTSV=NPRINT
      NPRINT=-5
C
      IF(.NOT.UNCON) LHNRUNC=LHNR
C     POINT LHNRUNC TO LHNR TO AVOID COPYING
C
C     NOTE THAT IN CASE OF NUMERIC HESSIANS ONE HAS TO TURN OFF
C     SYMMETRY.  HERE WE ONLY USE TOTALLY SYMMETRIC STRECHES, SO
C     SYMMETRY CAN BE USED IN ITS FULL GLORY!
C
      CALL TSECND(T0)
      IF(SOME) WRITE(IW,9000) NMODE,QRTOL,NVIB,NEED+NEEDP
C
C     ----- CALCULATE DERIVATIVE OF RELATIVISTIC INTEGRALS -----
C           BY NUMERICAL DIFFERENCE
C
      CALL VCLR(X(LDE),1,3*NAT)
      CALL DCOPY(3*NAT,C,1,X(LGE),1)
      IF(NVIB.EQ.1) THEN
C        SAVE THE RELATIVISTIC INTEGRALS AT THE INITIAL GEOMETRY
         CALL DAREAD(IDAF,IODA,X(LHNRUNC),LU2,NDARELB+6,0)
         CALL DAREAD(IDAF,IODA,X(LTAO2),LU2,NDARELB+4,0)
         CALL DAXPY(LU2,HALF,X(LTAO2),1,X(LHNRUNC),1)
         IF (UNCON) CALL CONTRAM(X(LHNR),X(LHNRUNC),LU2,0,1)
         CALL DAREAD(IDAF,IODA,X(LHR),L2,11,0)
         CALL VSUB(X(LHNR),1,X(LHR),1,X(LQRINT+L2),1,L2)
         CALL DSCAL(L2,-ONE/QRTOL,X(LQRINT+L2),1)
      ENDIF
C
C     PARALLEL RUNS ARE EXCRUTIATING PAIN IN THE MID-BODY FAT RESERVOIR.
C
      SVDSKW=DSKWRK
      SVGPAR=GOPARR
C
C     CREATE THE SHADOW DICTIONARY FILE. IT IS ABSOLUTELY NECESSARY
C     FOR PARALLEL RUNS DUE TO THE DAF MASTER-ONLY POLICY. HOWEVER,
C     IT IS TIME-SAVING TO DO THIS FOR SEQUENTIAL RUNS ALSO AS IT
C     AVOIDS REWRITING AND HAVING TO RECOMPUTE INTEGRALS AT THE
C     INITIAL GEOMETRY.
C     LENRA DOES NOT SEEM TO HAVE MUCH SIGNIFICANCE AT PRESENT
C
      IUDA=30
      LENRA=LU2
      DSKWRK=.TRUE.
      CALL RAOPEN(IUDA,X(LIUDA),0,NIODA,LENRA,NPRINT)
C
C        COPY RECORDS. IT IS AN *ERROR* TO COPY NONEXISTENT RECORDS.
C        SAVED ARE: CONTRACTED/UNCONTRACTED BASIS SETS AND THE
C        CONTRACTING MATRIX THAT CONNECTS THE TWO.
C        NONE OF THAT IS NEEDED FOR NO UNCONTRACTION.
C
      IF(UNCON) THEN
         CALL DAREAD(IDAF,IODA,    X(LBUF),NBUF1,NDARELB+15,0)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),NBUF1,NDARELB+15,0)
         CALL DAREAD(IDAF,IODA,    X(LBUF),NBUF2,NDARELB+16,1)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),NBUF2,NDARELB+16,1)
         CALL DAREAD(IDAF,IODA,    X(LBUF),NBUF1,NDARELB+17,0)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),NBUF1,NDARELB+17,0)
         CALL DAREAD(IDAF,IODA,    X(LBUF),NBUF2,NDARELB+18,1)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),NBUF2,NDARELB+18,1)
         CALL DAREAD(IDAF,IODA,    X(LBUF),NBUF3,NDARELB+19,0)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),NBUF3,NDARELB+19,0)
      ENDIF
      IF(KAMERAD.EQ.1.AND..NOT.UNCON) THEN
         CALL DAREAD(IDAF,IODA,    X(LBUF),L3,   44,0)
         CALL RAWRIT(IUDA,X(LIUDA),X(LBUF),L3,   44,0)
      ENDIF
      IDAFS=IDAF
      CALL ICOPY(NIODA,IODA,1,X(LJUDAS),1)
      JRECLN=IRECLN
      JRECST=IRECST
      CALL ICOPY(NIODA,IFILEN,1,X(LFILEN),1)
      IDAF=IUDA
      CALL ICOPY(NIODA,X(LIUDA),1,IODA,1)
C
C        SO, FROM NOW ON THE DICTIONARY FILE IS SECRETLY REPLACED.
C        DAOPEN/DAWRIT ARE IMPLANTED WITH A CALL TO RAREAD/RAWRIT
C
      NXT = IBTYP.EQ.1
      NEXT = -1
      MINE = -1
C
      DO 300 IMODE=1,NMODE
         SKIP=.FALSE.
         IF(GOPARR.AND.NXT) THEN
           MINE = MINE + 1
           IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
           IF(NEXT.NE.MINE) SKIP=.TRUE.
         END IF
         IF(GOPARR  .AND.  .NOT.NXT  .AND.  MOD(IMODE,NPROC).NE.ME)
     *        SKIP=.TRUE.
         IF(SKIP) THEN
           IF(SOME) WRITE(IW,9050) IMODE
           GOTO 300
         ENDIF
C
         IF(NVIB.EQ.2) THEN
            CALL VCLR(X(LDQRINT),1,L2)
         ELSE
            CALL DCOPY(L2,X(LQRINT+L2),1,X(LDQRINT),1)
         ENDIF
C
C        NO PARALLELISATION BEYOND THIS POINT (ONE NODE DOES TWO
C        OFFSETS).  ONEEI HAS PRIMITIVE INTEGRALS PARALLELISED.
C        THAT MUST BE STOPPED.
C
         GOPARR=.FALSE.
C
C     ------- MOVE IN THE +/- DIRECTIONS
C
         DO IVIB=1,NVIB
           DIFF=SSIGN(IVIB)*QRTOL
C
C          SHIFT EACH ATOM ALONG THE COMPUTED STRETCH
C
           CALL DAXPY(3*NAT,DIFF,X(LP+IMODE-1),NCCF,C,1)
C
           NPRSAVE = NPRINT
           NPRINT = -2183
           CALL ONEEI
           NPRINT = NPRSAVE
C
           CALL DAREAD(IDAF,IODA,X(LHNRUNC),LU2,NDARELB+6,0)
           CALL DAREAD(IDAF,IODA,X(LTAO2),LU2,NDARELB+4,0)
C
C          HNR=V+T
C
           CALL DAXPY(LU2,HALF,X(LTAO2),1,X(LHNRUNC),1)
           IF (UNCON) CALL CONTRAM(X(LHNR),X(LHNRUNC),LU2,0,1)
C
C          CONVERT NR INTEGRALS FROM UNC. BASES TO CONT. BASES
C
           CALL DAREAD(IDAF,IODA,X(LHR),L2,11,0)
C
C          VR=HNR-HR
C
           CALL VSUB(X(LHNR),1,X(LHR),1,X(LQRINT),1,L2)
C
C          DVR/DR
C
           CALL DAXPY(L2,ONE/(NVIB*DIFF),X(LQRINT),1,X(LDQRINT),1)
C
C          MOTO NI MODOSU (RETURN TO FORMER COORDINATES
C
           CALL DCOPY(3*NAT,X(LGE),1,C,1)
C
           IF(SOME) THEN
             CALL TSECND(T1)
             TCPU = T1 - T0
             T0 = T1
             WRITE(IW,9040) IMODE,IVIB,TCPU
           END IF
         ENDDO
C
         GOPARR=SVGPAR
C
C        CONTRIBUTION TO DE/DR:  D*DVR/DR. DVR/DR IS COMPUTED
C        NUMERICALLY (DERIVATIVE OF VR, THE RELATIVISTIC CONTRIBUTION
C        TO THE HAMILTONIAN) D IS DENSITY.
C
         IJ = 0
         SUM=0.0D+00
         DO I = 1, L1
           DO J = 1, I
             IJ = IJ + 1
C            DE(IXYZ,IATM)=DE(IXYZ,IATM)+DENAB(IJ)*X(LDQRINT-1+IJ)
             SUM=SUM+DENAB(IJ)*X(LDQRINT-1+IJ)
           END DO
         END DO
         X(LDE+IMODE-1)=SUM
         IF(DBG) WRITE(IW,9045) SUM*QRTOL,SUM*QRTOL/ETOT*1.0D+06
  300 CONTINUE
C
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
C
C    THERE IS NO NEED TO UPDATE DICTIONARY FILE, ONLY RESTORE ITS IODA
C    FROM NOW ON THE DICTIONARY FILE IS BACK TO NORMAL
C
      CALL RACLOS(IUDA,'DELETE')
      DSKWRK=SVDSKW
      IDAF=IDAFS
      CALL ICOPY(NIODA,X(LJUDAS),1,IODA,1)
      IRECLN=JRECLN
      IRECST=JRECST
      CALL ICOPY(NIODA,X(LFILEN),1,IFILEN,1)
C
C     CONVERT THE GRADIENT BACK TO CARTESIANS
C
      CALL MRARBR(X(LDE),1,1,NMODE,X(LP),NCCF,NAT*3,X(LDEC),1)
C
C     ADD THE RELATIVISTIC CONTRIBUTION TO THE TOTAL GRADIENT
C
      CALL DAXPY(3*NAT,1.0D+00,X(LDEC),1,DE,1)
      IF(DBG) THEN
         IF(GOPARR) CALL DDI_GSUMF(1502,X(LDEC),3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,*) 'CONTRIBUTION FROM THE RELATIVISTIC GRADIENT'
            CALL EGOUT(X(LDEC),NAT)
         ENDIF
      ENDIF
C
C     ----- FINALLY, RE-CALCULATE RELATIVISTIC INTEGRALS -----
C           THIS IS NEEDED FOR ONLY LAST OPTIMIZATION PROCESS.
C           PARALLEL RUNS DO NOT DESTROY INTEGRALS ON FILE 10 SO THEY
C           DO NOT HAVE TO BE RECOMPUTED.
C
C     IF(.NOT.GOPARR) THEN
C     IF(MASWRK) WRITE(IW,9100)
C9100 FORMAT(1X,'RECOMPUTING ONE-ELECTRON INTEGRALS FOR IVIB=0.',/)
C     CALL ONEEI
C     ENDIF
C
      CALL RETFM(NEED)
      CALL RETFM(NEEDP)
      NPRINT=NPRINTSV
C
      RETURN
 8000 FORMAT(//1X,62(1H!),/1X,'THE NUMERIC GRADIENT DOES NOT WORK WELL',
     *       ' WITH THE HONDO OPTION.',/1X,62(1H!),//)
 9000 FORMAT(/1X,71(1H-)/
     *        1X,'RELATIVISTIC 1E- GRADIENT WILL BE',
     *           ' COMPUTED BY NUMERICAL DIFFERENTIATION'/
     *           I4,' TOTALLY SYMMETRIC MODES FOUND,',
     *              ' EACH MUST BE DISPLACED'/
     *        1X,'COORDINATE OFFSET (BOHR):',F15.7,
     *           '  (',I1,' OFFSETS/MODE)'/
     *        1X,'MEMORY REQUIRED (WORDS):',I18/1X,71(1H-)/)
 9040 FORMAT(1X,'IMODE=',I3,' IVIB=',I1,
     *       ' ONE ELECTRON INTEGRALS TOOK ',F9.2,' SEC.')
 9045 FORMAT(9X,'ENERGY DIFFERENCE: ',F13.10,' (',F9.7,
     *          ' PPM OF THE TOTAL ENERGY)')
 9050 FORMAT(1X,'IMODE=',I3,' IS ASSIGNED TO ANOTHER NODE.')
      END
C*MODULE QREL     *DECK UNLINDEP
      SUBROUTINE UNLINDEP(N,M,V,LDV,T,LDT,WRK,TOL,NLI)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(LDV,M),WRK(N),T(LDT,M)
      PARAMETER (ONE=1.0D+00)
C
C     USE GRAM-SCHMIDT ORTHOGONALISATION TO PICK ONLY LINEARLY
C     INDEPENDENT (LI) VECTORS FROM -M- VECTORS -V- OF LENGTH -N-.
C     -T- WILL HOLD GRAM-SCHMIDT ORTHOGONALISED LI VECTORS AT THE END.
C     VECTORS IN -V- WILL HAVE THE LINEARLY DEPENDENT ONES SHIFTED
C     TO THE END, WITH -NLI- HOLDING THE NUMBER OF LINEARLY
C     INDEPENDENT VECTORS FOUND.  THE FIRST NLI VECTORS IN -V-
C     WILL BE ORTHONORMALIZED JUST BEFORE EXIT, AND MAY BE USED
C     TO SPAN THE "TOTALLY SYMMETRIC" VIBRATIONAL MODE SPACE.
C
C     TOL IS A CRITERION TO DECIDE LINEAR DEPENDENCY (VECTOR LENGTH).
C     INPUT: N,M,V,LDV,LDT,TOL
C     OUTPUT: V,T,NLI
C
      NLD=0
      NLI=0
      I=1
  100 CONTINUE
         IF(I.GT.M-NLD) GOTO 200
         CALL DCOPY(N,V(1,I),1,WRK,1)
         DO J=1,NLI
            CALL DAXPY(N,-DDOT(N,V(1,I),1,T(1,J),1),T(1,J),1,WRK,1)
         ENDDO
         Q=SQRT(DDOT(N,WRK,1,WRK,1))
         IF(ABS(Q).LT.TOL) THEN
            CALL DSWAP(N,V(1,I),1,V(1,M-NLD),1)
            NLD=NLD+1
         ELSE
            CALL DSCAL(N,ONE/Q,WRK,1)
            NLI=NLI+1
            CALL DCOPY(N,WRK,1,T(1,NLI),1)
            I=I+1
         END IF
         GOTO 100
  200 CONTINUE
C
C        SCHMIDT ORTHOGONALIZE LINEARLY INDEPENDENT VECTORS JUST FOUND
C
      DO I=1,NLI
         DO J=1,I-1
            DUM = DDOT(N,V(1,I),1,V(1,J),1)
            CALL DAXPY(N,-DUM,V(1,J),1,V(1,I),1)
         ENDDO
         DUM = ONE/SQRT(DDOT(N,V(1,I),1,V(1,I),1))
         CALL DSCAL(N,DUM,V(1,I),1)
      ENDDO
      RETURN
      END
C*MODULE QREL     *DECK DCLEAN
      SUBROUTINE  DCLEAN(N,DX,INCX,THR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DX(INCX,*)
C
C     REPLACE TINY NUMBERS IN DX BY ZEROS
C     TINY IS LESS THAN THR
C
      DO I=1,N
         IF(ABS(DX(1,I)).LT.THR) DX(1,I)=0.0D+00
      ENDDO
      RETURN
      END
C*MODULE QREL     *DECK QKIN2
      SUBROUTINE  QKIN2(XX,UU,TD,L0,L1,L2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XX(L1,L0),UU(L1,L0),TD(L0)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FMCOM / X(1)
      PARAMETER(NDQ=2,TWO=2.0D+00)
C
      CALL VALFM(LOADFM)
      LT=LOADFM+1
      LU=LT+L2*NDQ
      LTT=LU+L1*L1*NDQ
      LWRK=LTT+L2*NDQ
      LAST=LWRK+L1*NDQ
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
      CALL DAREAD(IDAF,IODA,X(LT),L2*NDQ,395,0)
      CALL DQCOPY(L1*L0,XX,1,X(LU),1)
      CALL QTFTRI(X(LTT),X(LT),X(LU),X(LWRK),L0,L1,L1)
      CALL DQDIAG(X(LTT),UU,TD,L1,L0,1)
C           WE NEED TO DIAGONALISE 2T, NOT T
      CALL DSCAL(L0,TWO,TD,1)
      CALL RETFM(NEED)
      RETURN
      END
