C
C $Id: pmix.F,v 1.2 1998/07/16 16:40:22 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE PMIX(ITER)
C
C     USES THE IMPROVED ITERATION SCHEME (IIS) TO ACCELERATE CONVERGENCE
C     ON THE DENSITY MATRIX.  FOR DETAILS SEE:  P. BADZIAG AND F. SOLMS,
C     COMPUT. CHEM., VOL 12, PP 233-236 (1988).  THE ITERATION COUNTER
C     CORRESPONDS TO THE SCF CYCLE, I.E., ITER=1 FOR THE FIRST ITERATION,
C     ITER=2 FOR THE SECOND, ETC.  THIS DIFFERS FROM BADZIAG AND SOLMS
C     CONVENTION WHERE ITER=0 CORRESPONDS TO THE FIRST CYCLE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
      LOGICAL MIXED
      SAVE A1,A2,XLOLD,XLNEW,MIXED
C
      IF(ITER.EQ.1)THEN
        A1 = 1.0D0
        A2 = 1.0D0
        XLNEW = 1.0D0
        MIXED = .FALSE.
      ENDIF
      IIMAX = IIMAT(NATOMS+1)-1
      IJMAX = IJMAT(IP1(NATOMS+1))-1
      IF(ITER.EQ.2)THEN
        DO 10 II=1,IIMAX
          DPII(II) = PDIAG(II) - PIIOLD(II)
 10     CONTINUE
      ENDIF
      IF(ITER.LT.3) RETURN
C
C     DO THE FIRST MIXING STEP IF WE'VE COMPLETED TWO ITERATIONS.
C
      IF(ITER.GE.3)THEN
        POLD = 1.0D0 - A2
        PNEW = A2
C-RDC c       write(0,*) ' first mix, POLD,PNEW=',POLD,PNEW
C-RDC c       write(8,*) ' first mix, POLD,PNEW=',POLD,PNEW
        DO 30 II=1,IIMAX
          PDIAG(II) = PNEW*PDIAG(II) + POLD*PIIOLD(II)
 30     CONTINUE
        DO 40 IJ=1,IJMAX
          PDIAT(IJ) = PNEW*PDIAT(IJ) + POLD*PIJOLD(IJ)
 40     CONTINUE
      ENDIF
C
C     GET NEW LAMBDA PARAMETER (XLNEW) AND FORM NEW DIFFERENCE MATRIX.
C
      XLOLD = XLNEW
      DD12 = 0.0D0
      DD11 = 0.0D0
      DO 70 I=1,NATOMS
        NORBSI = NATORB(IATNUM(I))
        IF(NORBSI.EQ.0) GO TO 70
        II = IIMAT(I)
        DO 60 IP=1,NORBSI
          DO 50 JP=1,IP
            DII = DPII(II)
            D11TMP = DII*DII
            DD11 = DD11 + D11TMP
            DPII(II) = PDIAG(II) - PIIOLD(II)
            D12TMP = DII*DPII(II)
            DD12 = DD12 + D12TMP
            IF(IP.EQ.JP)THEN
              DD11 = DD11 - 0.5D0*D11TMP
              DD12 = DD12 - 0.5D0*D12TMP
            ENDIF
            II = II + 1
 50       CONTINUE
 60     CONTINUE
 70   CONTINUE
      XLNEW = DD12/DD11
C-RDC c     write(0,*) ' XLNEW=',XLNEW
C-RDC c     write(8,*) ' XLNEW=',XLNEW
      XLTEST = ABS(XLNEW-XLOLD)/MAX(ABS(XLOLD),ABS(XLNEW))
      IF(.NOT.MIXED)THEN
        XLTEST = ABS(XLNEW-XLOLD)/MAX(ABS(XLOLD),ABS(XLNEW))
        IF(XLTEST.LE.0.1D0)THEN
          MIXED = .TRUE.
        ELSE
          RETURN
        ENDIF
      ENDIF
      DENOM = 1.0D0 - XLNEW
      IF(ABS(DENOM).LT.0.01D0) DENOM = SIGN(0.01D0,DENOM)
      A1 = 1.0D0/DENOM
      IF((A1.GT.0.0D0.AND.A1.LE.1.0D0))THEN
        PNEW = A1
        POLD = 1.0D0 - PNEW
C-RDC c       write(0,*) ' second mix, POLD,PNEW=',POLD,PNEW
C-RDC c       write(8,*) ' second mix, POLD,PNEW=',POLD,PNEW
        DO 90 II=1,IIMAX
          PDIAG(II) = PNEW*PDIAG(II) + POLD*PIIOLD(II)
          DPII(II) = PNEW*DPII(II)
 90     CONTINUE
        DO 100 IJ=1,IJMAX
          PDIAT(IJ) = PNEW*PDIAT(IJ) + POLD*PIJOLD(IJ)
 100    CONTINUE
        XLNEW = PNEW*XLNEW
      ENDIF
      A2 = A1*A2
      IF(A2.LE.0.0D0.OR.A2.GT.1.0D0)THEN
        DENOM = 1.0D0 + A2
        IF(ABS(DENOM).LT.0.01D0) DENOM = SIGN(0.01D0,DENOM)
        A2 = 2.0D0*A2/DENOM
      ENDIF
      IF(ABS(A2).GT.1.5D0) A2 = SIGN(1.5D0,A2)
      RETURN
      END