      SUBROUTINE PULAY(F,P,N,FPPF,FOCK,EMAT,LFOCK,NFOCK,MSIZE,START,PL)
      implicit double precision (A-H,O-Z)
      dimension F(*), P(*), FPPF(*), FOCK(*)
      include 'MOP6FILE'
      logical START
************************************************************************
*
*   PULAY USES DR. PETER PULAY'S METHOD FOR CONVERGENCE.
*         A MATHEMATICAL DESCRIPTION CAN BE FOUND IN
*         "P. PULAY, J. COMP. CHEM. 3, 556 (1982).
*
* ARGUMENTS:-
*         ON INPUT F      = FOCK MATRIX, PACKED, LOWER HALF TRIANGLE.
*                  P      = DENSITY MATRIX, PACKED, LOWER HALF TRIANGLE.
*                  N      = NUMBER OF ORBITALS.
*                  FPPF   = WORKSTORE OF SIZE MSIZE, CONTENTS WILL BE
*                           OVERWRITTEN.
*                  FOCK   =      "       "              "         "
*                  EMAT   = WORKSTORE OF AT LEAST 15**2 ELEMENTS.
*                  START  = logical, = TRUE TO START PULAY.
*                  PL     = UNDEFINED ELEMENT.
*      ON OUTPUT   F      = "BEST" FOCK MATRIX, = LINEAR COMBINATION
*                           OF KNOWN FOCK MATRICES.
*                  START  = FALSE
*                  PL     = MEASURE OF NON-SELF-CONSISTENCY
*                         = [F*P] = F*P - P*F.
*
************************************************************************
      common /KEYWRD/ KEYWRD
      common /NUMCAL/ NUMCAL
      dimension EMAT(20,20), EVEC(1000), COEFFS(20)
      character*241 KEYWRD
      logical  DEBUG
      data ICALCN/0/
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
         MAXLIM=6
         DEBUG=(INDEX(KEYWRD,'DEBUGPULAY') .NE.0)
      end if
      if (START) then
         LINEAR=(N*(N+1))/2
         MFOCK=MSIZE/LINEAR
         if (MFOCK.GT.MAXLIM) MFOCK=MAXLIM
         if (DEBUG)
     1      write(MIOB,'('' MAXIMUM SIZE:'',I5)') MFOCK
         NFOCK=1
         LFOCK=1
         START=.FALSE.
      else
         if (NFOCK.LT.MFOCK)  NFOCK=NFOCK+1
         if (LFOCK.NE.MFOCK) then
            LFOCK=LFOCK+1
         else
            LFOCK=1
         end if
      end if
      LBASE=(LFOCK-1)*LINEAR
*
*   FIRST, STORE FOCK MATRIX FOR FUTURE REFERENCE.
*
      do 10 I=1,LINEAR
   10 FOCK((I-1)*MFOCK+LFOCK)=F(I)
*
*   NOW FORM /FOCK*DENSITY-DENSITY*FOCK/, AND STORE THIS IN FPPF
*
      call MAMULT(P,F,FPPF(LBASE+1),N,0.D0)
      call MAMULT(F,P,FPPF(LBASE+1),N,-1.D0)
*
*   FPPF NOW CONTAINS THE RESULT OF FP - PF.
*
      NFOCK1=NFOCK+1
      do 20 I=1,NFOCK
         EMAT(NFOCK1,I)=-1.D0
         EMAT(I,NFOCK1)=-1.D0
         EMAT(LFOCK,I)=DOT(FPPF((I-1)*LINEAR+1),FPPF(LBASE+1),LINEAR)
   20 EMAT(I,LFOCK)=EMAT(LFOCK,I)
      PL=EMAT(LFOCK,LFOCK)/LINEAR
      EMAT(NFOCK1,NFOCK1)=0.D0
      CONST=1.D0/EMAT(LFOCK,LFOCK)
      do 30 I=1,NFOCK
         do 30 J=1,NFOCK
   30 EMAT(I,J)=EMAT(I,J)*CONST
      if (DEBUG) then
         write(MIOB,'('' EMAT'')')
         do 40 I=1,NFOCK1
   40    write(MIOB,'(6E13.6)')  (EMAT(J,I),J=1,NFOCK1)
      end if
      L=0
      do 50 I=1,NFOCK1
         do 50 J=1,NFOCK1
            L=L+1
   50 EVEC(L)=EMAT(I,J)
      CONST=1.D0/CONST
      do 60 I=1,NFOCK
         do 60 J=1,NFOCK
   60 EMAT(I,J)=EMAT(I,J)*CONST
*********************************************************************
*   THE MATRIX EMAT SHOULD HAVE FORM
*
*      |<E(1)*E(1)>  <E(1)*E(2)> ...   -1.0|
*      |<E(2)*E(1)>  <E(2)*E(2)> ...   -1.0|
*      |<E(3)*E(1)>  <E(3)*E(2)> ...   -1.0|
*      |<E(4)*E(1)>  <E(4)*E(2)> ...   -1.0|
*      |     .            .      ...     . |
*      |   -1.0         -1.0     ...    0. |
*
*   WHERE <E(I)*E(J)> IS THE SCALAR PRODUCT OF [F*P] FOR ITERATION I
*   TIMES [F*P] FOR ITERATION J.
*
*********************************************************************
      call OSINV(EVEC,NFOCK1,D)
      if (ABS(D).LT.1.D-6) then
         START=.TRUE.
         return
      end if
      if (NFOCK.LT.2) return
      IL=NFOCK*NFOCK1
      do 70 I=1,NFOCK
   70 COEFFS(I)=-EVEC(I+IL)
      if (DEBUG) then
         write(MIOB,'('' EVEC'')')
         write(MIOB,'(6F12.6)') (COEFFS(I),I=1,NFOCK)
         write(MIOB,'(''    LAGRANGIAN MULTIPLIER (ERROR) =''
     %                ,F13.6)')  EVEC(NFOCK1*NFOCK1)
      end if
      do 90 I=1,LINEAR
         SUM=0
         L=0
         II=(I-1)*MFOCK
         do 80 J=1,NFOCK
   80    SUM=SUM+COEFFS(J)*FOCK(J+II)
   90 F(I)=SUM
      return
      end
