subroutine mcweeny(nmcw)
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     This subroutine performs the McWeeny transformation
C     P' = 3P**2 - 2P**3
C
C     The old density matrix (stored in pdiag and pdiat) will
C     be overwritten with the transformed matrix. In the process
C     piiold and pijold will be used as storage space; piiold
C     and pijold will also be overwritten with the transformed
C     matrix.
C
C     The transformation is performed nmcw times.
C
C     Written by Arjan van der Vaart, 3/18/1999
C

C     to save on memory, store the rows in ff, and the columns in ww
      dimension b(4*maxatm,4), tmp(4,4)
      logical igtj, kgtj

C     we need the "true" density matrix, i.e. half the
C     closed shell density matrix
      npairs = ip1(natoms+1)-1
      iimax = iimat(natoms+1)-1
      ipmax = ip1(natoms+1)
      ijmax = ijmat(ipmax)-1
      do ii=1,iimax
         pdiag(ii) = 0.5*pdiag(ii)
      enddo
      do ij=1,ijmax
         pdiat(ij) = 0.5*pdiat(ij)
      enddo

      imcw = 0
 10   imcw = imcw + 1
        write(*,*)'mcweeny'
      if (imcw.gt.nmcw) goto 6000

C     To save memory, the tranformation is performed as:
C     P'(kat,jat) = Sum(i) [ P(kat,iat) . B(iat,jat) ]
C
C     with B(iat,jat) = 3 P(iat,jat) - 2 P**2(iat,jat)
C
C     This is also a very efficient way to calculate P';
C     the algorithm used is a fast block-algorithm which
C     minimizes the number of computations and retrieves/stores the
C     data sequentially.
C
      do 5000 kat=1,natoms

         katn = iatnum(kat)
         if (katn.eq.0) goto 5000
         norbsk = natorb(katn)

         do 4000 jat=1,kat

            jatn = iatnum(jat)
            if (jatn.eq.0) goto 4000
            norbsj = natorb(jatn)

            if (kat.ne.jat) then
               call ijfind(npairs,kat,jat,ijaddr)
               if (ijaddr.eq.0) goto 4000
C     .        kj will be used later on
               kj = ijmat(ijaddr)-1
            endif

            call getrow(npairs,jat,norbsj,ncolumn,ww)

C     .     start the loop calculating B(:,jat)
            ib0 = 0
            do 1000 iat=1,natoms
               iatn = iatnum(iat)
               if (iatn.eq.0) goto 1000
               norbsi = natorb(iatn)
               call getrow(npairs,iat,norbsi,nrow,ff)

               do i=1,4
                  do j=1,4
                     tmp(j,i) = 0.0
                  enddo
               enddo

C     .        calculate B(iat,jat)

               if (iat.gt.jat) then
                  igtj = .true.
                  i = iat
                  j = jat
               elseif (iat.eq.jat) then
                  iidiag = iimat(iat)-1
                  do i=1,norbsi
                     i0 = i*natoms4 - natoms4
                     do j=1,i
                        j0 = j*natoms4 - natoms4
                        do k=1,nrow
                           tmp(j,i) = tmp(j,i) + ff(k+i0)*ff(k+j0)
                        enddo
                        iidiag = iidiag + 1
                        tmp(j,i) = 3.0*pdiag(iidiag)-2.0*tmp(j,i)
                        tmp(i,j) = tmp(j,i)
                     enddo
                  enddo
                  goto 50
               else
                  igtj = .false.
                  i = jat
                  j = iat
               endif
               call ijfind(npairs,i,j,ijaddr)
               if (ijaddr.eq.0) goto 1000
               ijstart = ijmat(ijaddr)-1
               if (igtj) then
                  do j=1,norbsi
                     iy = ijstart+(j-1)*norbsj
                     j0 = j*natoms4 - natoms4
                     do i=1,norbsj
                        i0 = i*natoms4 - natoms4
                        do k=1,nrow
                           tmp(j,i) = tmp(j,i) + ff(k+j0)*ww(k+i0)
                        enddo
                        ix = iy+i
                        tmp(j,i) = 3.0*pdiat(ix)-2.0*tmp(j,i)
                     enddo
                  enddo
               else
                  do i=1,norbsj
                     ix = ijstart+(i-1)*norbsi
                     i0 = i*natoms4 - natoms4
                     do j=1,norbsi
                        j0 = j*natoms4 - natoms4
                        do k=1,nrow
                           tmp(j,i) = tmp(j,i) + ff(k+j0)*ww(k+i0)
                        enddo
                        iy = ix + j
                        tmp(j,i) = 3.0*pdiat(iy)-2.0*tmp(j,i)
                     enddo
                  enddo
               endif

 50            continue
               do j=1,norbsj
                  do i=1,norbsi
                     b(i+ib0,j) = tmp(i,j)
                  enddo
               enddo
               ib0 = ib0 + norbsi

 1000       continue

C     .     get P(kat,:)
            call getrow(npairs,kat,norbsk,nrow,ff)
            do i=1,4
               do j=1,4
                  tmp(j,i) = 0.0
               enddo
            enddo

C     .     calculate P(kat,iat)*B(iat,jat)
            do k=1,norbsk
               k0 = k*natoms4 - natoms4
               do j=1,norbsj
                  do l=1,nrow
                     tmp(j,k) = tmp(j,k) + b(l,j)*ff(l+k0)
                  enddo
               enddo
            enddo

C     .     store P'(kat,jat)
            if (kat.eq.jat) then
               iidiag = iimat(kat)-1
               do k=1,norbsk
                  do j=1,k
                     iidiag = iidiag + 1
                     piiold(iidiag) = tmp(j,k)
                  enddo
               enddo
            else
               do k=1,norbsk
                  do j=1,norbsj
                     kj = kj + 1
                     pijold(kj) = tmp(j,k)
                  enddo
               enddo
            endif

 3000       continue
 4000    continue
 5000 continue
         
C     repeat this nmcw times
      goto 10

C     now copy piiold and pijold back to pdiag and pdiat,
C     keeping in mind that we now need the closed shell
C     density matrix again.

 6000 do ii=1,iimax
         pdiag(ii) = 2.0*piiold(ii)
      enddo
      do ij=1,ijmax
         pdiat(ij) = 2.0*pijold(ij)
      enddo

      end

CC------------------------------------------------------CC

      subroutine getrow(npairs,iat,norbsi,nrow,row)
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     This algorithm retrieves the rows of the density matrix
C     for atom iat. The number of rows equals the number of
C     orbitals for atom iat.
C
C     Arjan van der Vaart, 3/18/1999

      dimension row(*)
      logical igtk

      nrow = 0
      do 1000 kat=1,natoms
         katn = iatnum(kat)
         if (katn.eq.0) goto 1000
         norbsk = natorb(katn)
         if (iat.gt.kat) then
            igtk = .true.
            i = iat
            j = kat
         elseif (iat.eq.kat) then
            iidiag = iimat(iat)-1
            n = nrow-natoms4
            do i=1,norbsi
               i0 = i*natoms4 + n
               in0 = i + n
               do j=1,i
                  iidiag = iidiag + 1
                  row(j+i0) = pdiag(iidiag)
                  row(in0+j*natoms4) = pdiag(iidiag)
               enddo
            enddo
            nrow = nrow + norbsi
            goto 1000
         else
            igtk = .false.
            i = kat
            j = iat
         endif
         call ijfind(npairs,i,j,ijaddr)
         if (ijaddr.eq.0) goto 1000
         ijstart = ijmat(ijaddr)-1
         if (igtk) then
            n = nrow-natoms4
            do j=1,norbsi
               iy = ijstart+(j-1)*norbsk
               j0 = j*natoms4 + n
               do i=1,norbsk
                  ix = iy+i
                  row(j0+i) = pdiat(ix)
               enddo
            enddo
            nrow = nrow + norbsk
         else
            n = nrow-natoms4
            do i=1,norbsk
               ix = ijstart+(i-1)*norbsi
               i0 = i+n
               do j=1,norbsi
                  iy = ix + j
                  row(i0+j*natoms4) = pdiat(iy)
               enddo
            enddo
            nrow = nrow + norbsk
         endif

 1000 continue

      end