      SUBROUTINE CNVG(PNEW, P, P1,NORBS, NITER, PL)
      implicit double precision (A-H,O-Z)
      include 'MOP6FILE'
      dimension P1(*), P(*), PNEW(*)
      logical EXTRAP
C***********************************************************************
C
C  CNVG IS A TWO-POINT INTERPOLATION ROUTINE FOR SPEEDING CONVERGENCE
C       OF THE DENSITY MATRIX.
C
C ON OUTPUT P      = NEW DENSITY MATRIX
C           P1     = DIAGONAL OF OLD DENSITY MATRIX
C           PL     = LARGEST DIFFERENCE BETWEEN OLD AND NEW DENSITY
C                    MATRIX DIAGONAL ELEMENTS
C***********************************************************************
      common/KEYWRD/ KEYWRD
      common /NUMCAL/ NUMCAL
      save RHFUHF
      character*241 KEYWRD
      data ICALCN/0/
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
         if (INDEX(KEYWRD,'UHF').NE.0) then
            RHFUHF=1.D0
         else
            RHFUHF=2.D0
         end if
      end if
      PL=0.0D00
      FACA=0.0D00
      DAMP=1.D10
      if (NITER.GT.3) DAMP=0.05D0
      FACB=0.0D00
      FAC=0.0D00
      II=MOD(NITER,3)
      EXTRAP=II.NE.0
      SUM1=0.D0
      K=0
      do 20 I=1,NORBS
         K=K+I
         A=PNEW(K)
         SUM1=SUM1+A
         SA=ABS(A-P(K))
         if (SA.GT.PL) PL=SA
         if (EXTRAP) go to 10
         FACA=FACA+SA**2
         FACB=FACB+(A-2.D00*P(K)+P1(I))**2
   10    P1(I)=P(K)
   20 P(K)=A
      if (FACB.LE.0.0D00) go to 30
      if (FACA.LT.(100.D00*FACB)) FAC=SQRT(FACA/FACB)
   30 IE=0
      SUM2=0.D0
      do 50 I=1,NORBS
         II=I-1
         do 40 J=1,II
            IE=IE+1
            A=PNEW(IE)
            P(IE)=A+FAC*(A-P(IE))
            PNEW(IE)=P(IE)
   40    continue
         IE=IE+1
         if (ABS(P(IE)-P1(I)) .GT. DAMP) then
            P(IE)=P1(I)+SIGN(DAMP,P(IE)-P1(I))
         else
            P(IE)=P(IE)+FAC*(P(IE)-P1(I))
         end if
         P(IE)=MIN(RHFUHF,MAX(P(IE),0.D0))
         SUM2=SUM2+P(IE)
   50 PNEW(IE)=P(IE)
C
C   RE-NORMALIZE IF ANY DENSITY MATRIX ELEMENTS HAVE BEEN TRUNCATED
C
      SUM0=SUM1
   60 if (SUM2.GT.1.D-3) then
         SUM=SUM1/SUM2
      else
         SUM=0.D0
      end if
      SUM1=SUM0
      if (SUM2.GT.1.D-3.AND.ABS(SUM-1.D0).GT.1.D-5) then
C#      write(MIOB,'(6F12.6)') (P((I*(I+1))/2),I=1,NORBS)
         SUM2=0.D0
         do 70 I=1,NORBS
            J=(I*(I+1))/2
C
C   ADD ON A SMALL NUMBER IN CASE AN OCCUPANCY IS EXACTLY ZERO
C
            P(J)=P(J)*SUM+1.D-20
            P(J)=MAX(P(J),0.D0)
C
C  SET UP RENORMALIZATION OVER PARTLY OCCUPIED M.O.'S ONLY.  FULL M.O.'S
C  CAN'T BE FILLED ANY MORE
C
            if (P(J).GT.RHFUHF) then
               P(J)=RHFUHF
               SUM1=SUM1-RHFUHF
            else
               SUM2=SUM2+P(J)
            end if
   70    PNEW(J)=P(J)
         go to 60
      end if
      return
      end
