      subroutine canonical(nmax,ncontractions,nocc,nconts,ncontp,ncontd,
     &                     nfuncatom,isave,fock,coeff,eigenvalue,
     &                     fockmat,wkmat1,wkmat2,veci,vecj)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension nfuncatom(*),isave(*)

      dimension fock(*),coeff(nmax,*),eigenvalue(*),
     &          fockmat(ncontractions,*),wkmat1(ncontractions,*),
     &          wkmat2(ncontractions,*),veci(*),vecj(*)

      data loopmax/100/

      data zero,half,one,two/0.0d0,0.5d0,1.0d0,2.0d0/
      data tolerance,big,ev/0.000000000001d0,10000.0d0,27.212d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      if(myid.eq.0) rewind 44

      do 1001 i=1,ncontractions
 1001 if(myid.eq.0) call fastrd(44,wkmat2(1,i),ncontractions)

      call mpi_bcast(wkmat2,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      do 1002 j=1,ncontractions
      do 1002 i=1,ncontractions
      ij=((max(i,j)*(max(i,j)-1))/2)+min(i,j)
 1002 fockmat(i,j)=fock(ij)

      do 1003 j=1,ncontractions
      do 1004 i=1,ncontractions
 1004 wkmat1(i,j)=zero
      do 1003 k=1,ncontractions
      do 1003 i=1,ncontractions
 1003 wkmat1(i,j)=wkmat1(i,j)+fockmat(i,k)*wkmat2(k,j)

      do 1005 j=1,ncontractions
      do 1006 i=1,ncontractions
 1006 fockmat(i,j)=zero
      do 1005 i=1,ncontractions
      do 1005 k=1,ncontractions
 1005 fockmat(i,j)=fockmat(i,j)+wkmat2(k,i)*wkmat1(k,j)

      do 1007 j=1,ncontractions
      do 1008 i=1,ncontractions
 1008 wkmat1(i,j)=zero
 1007 wkmat1(j,j)=one

      do 1009 loop=1,loopmax

      f_ij_max=zero

      do 1010 i=1,ncontractions
      do 1010 j=1,i

      if(i.eq.j) goto 1010

      if(abs(fockmat(i,j)).gt.f_ij_max) f_ij_max=abs(fockmat(i,j))

      if(abs(fockmat(i,j)).gt.tolerance) then

                    theta=(fockmat(j,j)-fockmat(i,i))/(two*fockmat(i,j))

                                         else

                    goto 1010

                                         endif

      if(theta.lt.zero) then
                              tan_phi=theta+sqrt(theta*theta+1)
                        else
                              tan_phi=theta-sqrt(theta*theta+1)
                        endif

      if(abs(theta).gt.big) tan_phi=-(half/theta)

      phi=atan(tan_phi)

      c=cos(phi)
      s=sin(phi)
      tau=tan(phi*half)

      f_ii_save=fockmat(i,i)
      f_jj_save=fockmat(j,j)
      f_ij_save=fockmat(i,j)

      do 1011 ii=1,ncontractions
 1011 veci(ii)=fockmat(ii,i)

      do 1012 jj=1,ncontractions
 1012 vecj(jj)=fockmat(jj,j)

      do 1013 ii=1,ncontractions
 1013 fockmat(ii,i)=c*veci(ii)+s*vecj(ii)

      do 1014 jj=1,ncontractions
 1014 fockmat(jj,j)=c*vecj(jj)-s*veci(jj)

      do 1015 ii=1,ncontractions
 1015 fockmat(i,ii)=fockmat(ii,i)

      do 1016 jj=1,ncontractions
 1016 fockmat(j,jj)=fockmat(jj,j)

      fockmat(i,i)=c*c*f_ii_save+two*c*s*f_ij_save+s*s*f_jj_save
      fockmat(j,j)=c*c*f_jj_save-two*c*s*f_ij_save+s*s*f_ii_save

      fockmat(i,j)=zero
      fockmat(j,i)=zero

      do 1017 ii=1,ncontractions
 1017 veci(ii)=wkmat1(ii,i)

      do 1018 jj=1,ncontractions
 1018 vecj(jj)=wkmat1(jj,j)

      do 1019 ii=1,ncontractions
 1019 wkmat1(ii,i)=c*veci(ii)+s*vecj(ii)

      do 1020 jj=1,ncontractions
 1020 wkmat1(jj,j)=c*vecj(jj)-s*veci(jj)

 1010 continue

      if(f_ij_max.lt.tolerance) goto 1021

      if((myid.eq.0).and.(loop.eq.loopmax)) write(6,1022) f_ij_max
 1022 format(/,' orthogonalization routine did not fully converge,',/,
     &         ' maximum value = ', f21.16,/)

      if(loop.eq.loopmax) stop

 1009 continue

 1021 continue

      do 1023 i=1,ncontractions
 1023 isave(i)=i

 1024 iswitch=0

      do 1025 i=2,ncontractions
      j=isave(i)
      k=isave(i-1)
      if(fockmat(j,j).lt.fockmat(k,k)) then
                                             iswitch=1
                                             isave(i)=k
                                             isave(i-1)=j
                                       endif
 1025 continue

      if(iswitch.eq.1) goto 1024

      do 1026 i=1,ncontractions
 1026 eigenvalue(i)=fockmat(isave(i),isave(i))

      do 1027 j=1,ncontractions
      do 1028 i=1,ncontractions
 1028 coeff(i,j)=zero
      do 1027 k=1,ncontractions
      do 1027 i=1,ncontractions
 1027 coeff(i,j)=coeff(i,j)+wkmat2(i,k)*wkmat1(k,isave(j))

      do 1029 i=1,ncontractions

      if((myid.eq.0).and.(i.le.nocc)) write(6,1030) i,ev*eigenvalue(i)
 1030 format(/,'     orbital # ',i5,
     &         ',     energy = ',f15.6,' eV    (occupied)',/)

      if((myid.eq.0).and.(i.gt.nocc)) write(6,1031) i,ev*eigenvalue(i)
 1031 format(/,'     orbital # ',i5,
     &         ',     energy = ',f15.6,' eV    (virtual)',/)

      ii=0

      do 1032 j=1,nconts
      if(myid.eq.0) write(6,1033) ii+1,nfuncatom(ii+1),coeff(ii+1,i)
 1033 format('         a.o. #',i5,',     center',i4,',    s  ',f15.8)
 1032 ii=ii+1

      do 1034 j=1,ncontp
      if(myid.eq.0) write(6,1035) ii+1,nfuncatom(ii+1),coeff(ii+1,i)
 1035 format('         a.o. #',i5,',     center',i4,',    px ',f15.8)
      if(myid.eq.0) write(6,1036) ii+2,nfuncatom(ii+2),coeff(ii+2,i)
 1036 format('         a.o. #',i5,',     center',i4,',    py ',f15.8)
      if(myid.eq.0) write(6,1037) ii+3,nfuncatom(ii+3),coeff(ii+3,i)
 1037 format('         a.o. #',i5,',     center',i4,',    pz ',f15.8)
 1034 ii=ii+3

      do 1038 j=1,ncontd
      if(myid.eq.0) write(6,1039) ii+1,nfuncatom(ii+1),coeff(ii+1,i)
 1039 format('         a.o. #',i5,',     center',i4,',    dxx',f15.8)
      if(myid.eq.0) write(6,1040) ii+2,nfuncatom(ii+2),coeff(ii+2,i)
 1040 format('         a.o. #',i5,',     center',i4,',    dxy',f15.8)
      if(myid.eq.0) write(6,1041) ii+3,nfuncatom(ii+3),coeff(ii+3,i)
 1041 format('         a.o. #',i5,',     center',i4,',    dxz',f15.8)
      if(myid.eq.0) write(6,1042) ii+4,nfuncatom(ii+4),coeff(ii+4,i)
 1042 format('         a.o. #',i5,',     center',i4,',    dyy',f15.8)
      if(myid.eq.0) write(6,1043) ii+5,nfuncatom(ii+5),coeff(ii+5,i)
 1043 format('         a.o. #',i5,',     center',i4,',    dyz',f15.8)
      if(myid.eq.0) write(6,1044) ii+6,nfuncatom(ii+6),coeff(ii+6,i)
 1044 format('         a.o. #',i5,',     center',i4,',    dzz',f15.8)
 1038 ii=ii+6

 1029 continue

      return
      end
