      SUBROUTINE INTERP(N,NP,NQ,MODE,E,FP,CP,VEC,FOCK,P,H,VECL)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension FP(MPACK), CP(N,N)
      dimension VEC(N,N), FOCK(N,N),
     1          P(N,N), H(N*N), VECL(N*N)
**********************************************************************
*
* INTERP: AN INTERPOLATION PROCEDURE FOR FORCING SCF CONVERGANCE
*         ORIGINAL THEORY AND FORTRAN WRITTEN BY R.N. CAMP AND
*         H.F. KING, J. CHEM. PHYS. 75, 268 (1981)
**********************************************************************
*
* ON INPUT N     = NUMBER OF ORBITALS
*          NP    = NUMBER OF FILLED LEVELS
*          NQ    = NUMBER OF EMPTY LEVELS
*          MODE  = 1, do NOT RESET.
*          E     = ENERGY
*          FP    = FOCK MATRIX, AS LOWER HALF TRIANGLE, PACKED
*          CP    = EIGENVECTORS OF FOCK MATRIX OF ITERATION -1
*                  AS PACKED ARRAY OF N*N COEFFICIENTS
*
* ON OUTPUT CP   = BEST GUESSED SET OF EIGENVECTORS
*           MODE = 2 OR 3 - USED BY callING PROGRAM
**********************************************************************
      dimension THETA(MAXORB), IA(MAXORB)
      common /KEYWRD/ KEYWRD
      common /NUMCAL/ NUMCAL
      common/FIT/NPNTS,IDUM2,XLOW,XHIGH,XMIN,EMIN,DEMIN,X(12),F(12),
     1 DF(12)
      logical DEBUG, DEBUG1
      character*241 KEYWRD
      save ZERO, FF, RADMAX, ICALCN, DEBUG1, DEBUG, IA
      data ICALCN/0/
      data ZERO,FF,RADMAX/0.0D0,0.9D0,1.5708D0/
      if (ICALCN.NE.NUMCAL) then
         DEBUG=(INDEX(KEYWRD,'INTERP').NE.0)
         DEBUG1=(INDEX(KEYWRD,'DEBUG').NE.0.AND.DEBUG)
         ICALCN=NUMCAL
         do 10 I=1,MAXORB
   10    IA(I)=(I*I-I)/2
      end if
C
C     RADMAX=MAXIMUM ROTATION ANGLE (RADIANS).  1.5708 = 90 DEGREES.
C         FF=FACTOR FOR CONVERGENCE TEST FOR 1D SEARCH.
C
      MINPQ=MIN0(NP,NQ)
      MAXPQ=MAX0(NP,NQ)
      NP1=NP+1
      NP2=MAX0(1,NP/2)
      if (MODE.EQ.2) go to 110
C
C     (MODE=1 OR 3 ENTRY)
C     TRANSFORM FOCK MATRIX TO CURRENT MO BASIS.
C     ONLY THE OFF DIAGONAL OCC-VIRT BLOCK IS COMPUTED.
C     STORE IN FOCK ARRAY
C
      II=0
      do 50 I=1,N
         I1=I+1
         do 40 J=1,NQ
            DUM=ZERO
            do 20 K=1,I
   20       DUM=DUM+FP(II+K)*CP(K,J+NP)
            if (I.EQ.N) go to 40
            IK=II+I+I
            do 30 K=I1,N
               DUM=DUM+FP(IK)*CP(K,J+NP)
   30       IK=IK+K
   40    P(I,J)=DUM
   50 II=II+I
      do 80 I=1,NP
         do 70 J=1,NQ
            DUM=ZERO
            do 60 K=1,N
   60       DUM=DUM+CP(K,I)*P(K,J)
            FOCK(I,J)=DUM
   70    continue
   80 continue
      if (MODE.EQ.3) go to 100
C
C     CURRENT POINT BECOMES OLD POINT (MODE=1 ENTRY)
C
      do 90 I=1,N
         do 90 J=1,N
   90 VEC(I,J)=CP(I,J)
      EOLD=E
      XOLD=1.0
      MODE=2
      return
C
C     (MODE=3 ENTRY)
C     FOCK CORRESPONDS TO CURRENT POINT IN CORRESPONDING REPRESENTATION.
C     VEC DOES NOT HOLD CURRENT VECTORS. VEC SET IN LAST MODE=2 ENTRY.
C
  100 NPNTS=NPNTS+1
      if (DEBUG)
     %  write(MIOB,'(''   INTERPOLATED ENERGY:'',F13.6)')  E*23.061D0
      IPOINT=NPNTS
      go to 500
C
C    (MODE=2 ENTRY) CALCULATE THETA, AND U, V, W MATRICES.
C                   U ROTATES CURRENT INTO OLD MO.
C                   V ROTATES CURRENT INTO CORRESPONDING CURRENT MO.
C                   W ROTATES OLD INTO CORRESPONDING OLD MO.
C
  110 J1=1
      do 140 I=1,N
         if (I.EQ.NP1) J1=NP1
         do 130 J=J1,N
            P(I,J)=ZERO
            do 120 K=1,N
  120       P(I,J)=P(I,J)+CP(K,I)*VEC(K,J)
  130    continue
  140 continue
C
C     U = CP(DAGGER)*VEC IS NOW IN P ARRAY.
C     VEC IS NOW AVAILABLE FOR TEMPORARY STORAGE.
C
      IJ=0
      do 170 I=1,NP
         do 160 J=1,I
            IJ=IJ+1
            H(IJ)=0.D0
            do 150 K=NP1,N
  150       H(IJ)=H(IJ)+P(I,K)*P(J,K)
  160    continue
  170 continue
      call HQRII(H,NP,NP,THETA,VECL)
      do 180 I=NP,1,-1
         IL=I*NP-NP
         do 180 J=NP,1,-1
  180 VEC(J,I)=VECL(J+IL)
      do 200 I=1,NP2
         DUM=THETA(NP1-I)
         THETA(NP1-I)=THETA(I)
         THETA(I)=DUM
         do 190 J=1,NP
            DUM=VEC(J,NP1-I)
            VEC(J,NP1-I)=VEC(J,I)
  190    VEC(J,I)=DUM
  200 continue
      do 210 I=1,MINPQ
         THETA(I)=MAX(THETA(I),ZERO)
         THETA(I)=MIN(THETA(I),1.D0)
  210 THETA(I)=ASIN(SQRT(THETA(I)))
C
C     THETA MATRIX HAS NOW BEEN CALCULATED, ALSO UNITARY VP MATRIX
C     HAS BEEN CALCULATED AND STORED IN FIRST NP COLUMNS OF VEC MATRIX.
C     NOW COMPUTE WQ
C
      do 240 I=1,NQ
         do 230 J=1,MINPQ
            VEC(I,NP+J)=ZERO
            do 220 K=1,NP
  220       VEC(I,NP+J)=VEC(I,NP+J)+P(K,NP+I)*VEC(K,J)
  230    continue
  240 continue
      call SCHMIT(VEC(1,NP1),NQ,N)
C
C     UNITARY WQ MATRIX NOW IN LAST NQ COLUMNS OF VEC MATRIX.
C     TRANSPOSE NP BY NP BLOCK OF U STORED IN P
C
      do 260 I=1,NP
         do 250 J=1,I
            DUM=P(I,J)
            P(I,J)=P(J,I)
  250    P(J,I)=DUM
  260 continue
C
C     CALCULATE WP MATRIX AND HOLD IN FIRST NP COLUMNS OF P
C
      do 300 I=1,NP
         do 270 K=1,NP
  270    H(K)=P(I,K)
         do 290 J=1,NP
            P(I,J)=ZERO
            do 280 K=1,NP
  280       P(I,J)=P(I,J)+H(K)*VEC(K,J)
  290    continue
  300 continue
      call SCHMIB(P,NP,N)
C
C     CALCULATE VQ MATRIX AND HOLD IN LAST NQ COLUMNS OF P MATRIX.
C
      do 340 I=1,NQ
         do 310 K=1,NQ
  310    H(K)=P(NP+I,NP+K)
         do 330 J=NP1,N
            P(I,J)=ZERO
            do 320 K=1,NQ
  320       P(I,J)=P(I,J)+H(K)*VEC(K,J)
  330    continue
  340 continue
      call SCHMIB(P(1,NP1),NQ,N)
C
C     CALCULATE (DE/DX) AT OLD POINT
C
      DEDX=ZERO
      do 370 I=1,NP
         do 360 J=1,NQ
            DUM=ZERO
            do 350 K=1,MINPQ
  350       DUM=DUM+THETA(K)*P(I,K)*VEC(J,NP+K)
  360    DEDX=DEDX+DUM*FOCK(I,J)
  370 continue
C
C     STORE OLD POINT INFORMATION FOR SPLINE FIT
C
      DEOLD=-4.0D0*DEDX
      X(2)=XOLD
      F(2)=EOLD
      DF(2)=DEOLD
C
C     MOVE VP OUT OF VEC ARRAY INTO FIRST NP COLUMNS OF P MATRIX.
C
      do 380 I=1,NP
         do 380 J=1,NP
  380 P(I,J)=VEC(I,J)
      K1=0
      K2=NP
      do 410 J=1,N
         if (J.EQ.NP1) K1=NP
         if (J.EQ.NP1) K2=NQ
         do 400 I=1,N
            DUM=ZERO
            do 390 K=1,K2
  390       DUM=DUM+CP(I,K1+K)*P(K,J)
  400    VEC(I,J)=DUM
  410 continue
C
C     CORRESPONDING CURRENT MO VECTORS NOW HELD IN VEC.
C     COMPUTE VEC(DAGGER)*FP*VEC
C     STORE OFF-DIAGONAL BLOCK IN FOCK ARRAY.
C
  420 II=0
      do 460 I=1,N
         I1=I+1
         do 450 J=1,NQ
            DUM=ZERO
            do 430 K=1,I
  430       DUM=DUM+FP(II+K)*VEC(K,J+NP)
            if (I.EQ.N) go to 450
            IK=II+I+I
            do 440 K=I1,N
               DUM=DUM+FP(IK)*VEC(K,J+NP)
  440       IK=IK+K
  450    P(I,J)=DUM
  460 II=II+I
      do 490 I=1,NP
         do 480 J=1,NQ
            DUM=ZERO
            do 470 K=1,N
  470       DUM=DUM+VEC(K,I)*P(K,J)
            FOCK(I,J)=DUM
  480    continue
  490 continue
C
C     SET LIMITS ON RANGE OF 1-D SEARCH
C
      NPNTS=2
      IPOINT=1
      XNOW=ZERO
      XHIGH=RADMAX/THETA(1)
      XLOW=-0.5D0*XHIGH
C
C     CALCULATE (DE/DX) AT CURRENT POINT AND
C     STORE INFORMATION FOR SPLINE FIT
C     ***** JUMP POINT FOR MODE=3 ENTRY *****
C
  500 DEDX=ZERO
      do 510 K=1,MINPQ
  510 DEDX=DEDX+THETA(K)*FOCK(K,K)
      DENOW=-4.0D0*DEDX
      ENOW=E
      if (IPOINT.LE.12) go to 530
  520 format(//,'EXCESSIVE DATA PNTS FOR SPLINE.',/
     1,'IPOINT =',I3,'MAXIMUM IS 12.')
C
C     PERFORM 1-D SEARCH AND DETERMINE EXIT MODE.
C
  530 X(IPOINT)=XNOW
      F(IPOINT)=ENOW
      DF(IPOINT)=DENOW
      call SPLINE
      if ((EOLD-ENOW).GT.FF*(EOLD-EMIN).OR.IPOINT.GT.10) go to 560
C
C     (MODE=3 EXIT) RECOMPUTE CP VECTORS AT PREDICTED MINIMUM.
C
      XNOW=XMIN
      do 550 K=1,MINPQ
         CK=COS(XNOW*THETA(K))
         SK=SIN(XNOW*THETA(K))
         if (DEBUG)
     %      write(MIOB,'('' ROTATION ANGLE:'',F12.4)') SK*57.29578
         do 540 I=1,N
            CP(I,K)   =CK*VEC(I,K)-SK*VEC(I,NP+K)
  540    CP(I,NP+K)=SK*VEC(I,K)+CK*VEC(I,NP+K)
  550 continue
      MODE=3
      return
C
C     (MODE=2 EXIT) CURRENT VECTORS GIVE SATISFACTORY ENERGY IMPROVEMENT
C     CURRENT POINT BECOMES OLD POINT FOR THE NEXT 1-D SEARCH.
C
  560 if (MODE.EQ.2) go to 580
      do 570 I=1,N
         do 570 J=1,N
  570 VEC(I,J)=CP(I,J)
      MODE=2
  580 ROLD=XOLD*THETA(1)*57.29578D0
      RNOW=XNOW*THETA(1)*57.29578D0
      RMIN=XMIN*THETA(1)*57.29578D0
      if (DEBUG) write(MIOB,600) XOLD,EOLD*23.061,DEOLD,ROLD
     1,                          XNOW,ENOW*23.061D0,DENOW,RNOW
     2,                          XMIN,EMIN*23.061D0,DEMIN,RMIN
      EOLD=ENOW
      if (NPNTS.LE.200) return
      write(MIOB,610)
      do 590 K=1,NPNTS
  590 write(MIOB,620) K,X(K),F(K),DF(K)
      write(MIOB,630)
      return
  600 format(
     1/14X,3H X ,10X,6H F(X) ,9X,7H DF/DX ,21H   ROTATION (DEGREES),
     2/10H      OLD ,F10.5,3F15.10,
     3/10H  CURRENT ,F10.5,3F15.10,
     4/10H PREDICTED,F10.5,3F15.10/)
  610 format(3H  K,10H     X(K) ,15H       F(K)    ,10H     DF(K))
  620 format(I3,F10.5,2F15.10)
  630 format(10X)
      end
