      subroutine cdpops(level,n,natoms,nmaxatoms,nconts,ncontp,ncontd,
     &                  nbsystems,nshels,nshelp,nsheld,icfunc,
     &                  nfuncatom,isys,overlap,dmta,dmtb,cdnorm,
     &                  pmat,smat,psmat,population)

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(*),
     &          isys(nmaxatoms,*)

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

      dimension cdnorm(*),pmat(n,*),smat(n,*),psmat(n,*),population(*)

      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 psmat(i,j)=zero
      do 1003 k=1,n
      do 1003 i=1,n
 1003 psmat(i,j)=psmat(i,j)+pmat(i,k)*smat(k,j)

      ii=0

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

      do 1006 i=1,ncontp
      do 1007 k=1,3
 1007 nfuncatom(ii+k)=icfunc(nshelp(i))
 1006 ii=ii+3

      do 1008 i=1,ncontd
      do 1009 k=1,6
 1009 nfuncatom(ii+k)=icfunc(nsheld(i))
 1008 ii=ii+6

      do 1010 i=1,natoms
      population(i)=zero
      do 1010 j=1,n
 1010 if(nfuncatom(j).eq.i) population(i)=population(i)+psmat(j,j)

      do 1011 i=1,nbsystems
 1011 cdnorm(i)=zero

      do 1012 i=1,nbsystems
      do 1012 j=1,isys(1,i)
      do 1012 k=1,level
 1012 cdnorm(i)=cdnorm(i)+population((k-1)*(natoms/level)+isys(j+1,i))

      return
      end
