      subroutine mulliken(n,natoms,nconts,ncontp,ncontd,nshels,nshelp,
     &                    nsheld,icfunc,nfuncatom,natomtype,overlap,
     &                    dmta,dmtb,pmat,smat,pstmat,pssmat,pop)

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)

      dimension nshels(*),nshelp(*),nsheld(*),nfuncatom(*),icfunc(*),
     &          natomtype(*)

      dimension overlap(*),dmta(*),dmtb(*)

      dimension pmat(n,*),smat(n,*),pstmat(n,*),pssmat(n,*),pop(*)

      data zero,half,one/0.0d0,0.5d0,1.0d0/

      ij=0

      do 1001 i=1,n
      do 1001 j=1,i
      ij=ij+1
      factor=half
      if(i.eq.j) factor=one
      pmat(i,j)=(dmta(ij)+dmtb(ij))*factor
 1001 pmat(j,i)=(dmta(ij)+dmtb(ij))*factor

      ij=0

      do 1002 i=1,n
      do 1002 j=1,i
      ij=ij+1
      smat(i,j)=overlap(ij)
 1002 smat(j,i)=overlap(ij)

      do 1003 j=1,n
      do 1004 i=1,n
 1004 pstmat(i,j)=zero
      do 1003 k=1,n
      do 1003 i=1,n
 1003 pstmat(i,j)=pstmat(i,j)+pmat(i,k)*smat(k,j)

      do 1005 i=1,n
      do 1005 j=1,i
      ij=ij+1
      factor=half
      if(i.eq.j) factor=one
      pmat(i,j)=(dmta(ij)-dmtb(ij))*factor
 1005 pmat(j,i)=(dmta(ij)-dmtb(ij))*factor

      do 1006 j=1,n
      do 1007 i=1,n
 1007 pssmat(i,j)=zero
      do 1006 k=1,n
      do 1006 i=1,n
 1006 pssmat(i,j)=pssmat(i,j)+pmat(i,k)*smat(k,j)

      ii=0

      do 1008 i=1,nconts
      nfuncatom(ii+1)=icfunc(nshels(i))
 1008 ii=ii+1

      do 1009 i=1,ncontp
      do 1010 k=1,3
 1010 nfuncatom(ii+k)=icfunc(nshelp(i))
 1009 ii=ii+3

      do 1011 i=1,ncontd
      do 1012 k=1,6
 1012 nfuncatom(ii+k)=icfunc(nsheld(i))
 1011 ii=ii+6

      write(6,1013)
 1013 format(//,' mulliken population analysis :',/)

      do 1014 i=1,natoms

      pop(i)=zero

      do 1015 j=1,n
 1015 if(nfuncatom(j).eq.i) pop(i)=pop(i)+pstmat(j,j)

      write(6,1016) i,real64(natomtype(i))-pop(i)
 1016 format(5x,' atom # ',i4,'       charge = ',f8.4)

 1014 continue

      write(6,1017)
 1017 format(//,' mayer bond order analysis :',/)

      do 1018 i=1,natoms
      do 1018 j=i,natoms

      if(i.eq.j) goto 1018

      bond=zero

      do 1019 k=1,n
      do 1019 l=1,n
 1019 if((nfuncatom(k).eq.i).and.
     &   (nfuncatom(l).eq.j)) bond=bond+pstmat(k,l)*pstmat(l,k)
     &                                 -pssmat(k,l)*pssmat(l,k)

      write(6,1020) i,j,bond
 1020 format(5x,' atoms # ',i4,' and ',i4,'       bond order = ',f8.4)

 1018 continue

      write(6,1021)
 1021 format(//,' mayer valence indices :',/)

      do 1022 i=1,natoms

      valence=pop(i)+pop(i)

      do 1023 j=1,n
      do 1023 k=1,n
 1023 if((nfuncatom(j).eq.i).and.
     &   (nfuncatom(k).eq.i)) valence=valence-pstmat(j,k)*pstmat(k,j)

      write(6,1024) i,valence
 1024 format(5x,' atom # ',i4,'       valence = ',f8.4)

 1022 continue

      return
      end
