*********************************************************************************
*     TOP_POP: purpose integrate density over elf basins
*     and calculate fluctuations of the basin populations
*     The wavefunction is read on a file generated by Gaussian92/94
*     with the option out=wfn
*     The grid point assignment is read on a file basin.sbf calculated by 
*     TOP_BAS
*
*     The input can be made either from the keybord or from an input file
*     Input description
*     1. wfn file name (a40)
*     2. Threshold for deleting negligible contributions to overlaps
*        (recommended 10-8)
*
*     References
*     A. Becke and K. E. Edgecombe, J. Chem. Phys., 92, 5397-5404 (1990)
*     B. Silvi and A. Savin, Nature, 371, 683-686 (1994)
*     A. Savin, B. Silvi and F. Colonna, Can. J. Chem., 74, 1088-1096 (1996)
*
*     please report bugs or problems via electronic mail at:
*     silvi@lct.jussieu.fr
*
********************************************************************************
      program top_pop
      use confs_module
      use grid_module
      use orbitals_module
      use popul_module
      use eval_module
      use basins_module
      implicit none
      double precision, dimension (:), allocatable :: volume,vol_at
      double precision, dimension (:), allocatable :: stdev,rho_bas,cc 
      double precision, dimension (:), allocatable :: rho_at,taa,tbb 
      double precision, dimension (:), allocatable :: tab,tbas,t_i,p_i 
      double precision, dimension (:), allocatable :: popaa,popbb,popab
      double precision, dimension (:), allocatable :: pop_orb,popa,popb
      double precision, dimension (:), allocatable :: ta,tb,sbas
      double precision, dimension (:,:), allocatable :: bij,yij,sij,tij
      double precision, dimension (:,:), allocatable :: aaij,bbij,abij
      double precision, dimension (:,:), allocatable :: stbas,cova,covb
      double precision, dimension (:,:), allocatable :: covab,p_ij,t_ij
      double precision, dimension (:,:), allocatable :: at_contrib
      double precision, dimension (:,:,:), allocatable :: stij
      integer*2, dimension (:), allocatable :: ind_at,temp_label
      integer*2 kind
      integer, dimension (:), allocatable :: induhfa,induhfb,indi
      integer, dimension (:), allocatable :: indice_x,indice_y,indice_z
      integer, dimension (:), allocatable :: indice_ij
      integer, dimension (:), allocatable :: elf_ind,aim_ind
      double precision rho_tol,tota,totb,arga,argb,prod_mo 
      double precision erf0(0:10000),erf1(0:10000),pi,arg,const
      double precision dx, dy, dz, derf
      double precision weight,units,total,rho_box
      double precision h_alpha,h_beta
      double precision pmatrix
      double precision sum,sum2
      double precision totals,total_charge
      double precision amult
      double precision tol
      double precision totala,totalb, totalaa,totalbb,totalab
      real t(2),spent,elapsed
      integer k,npt,nbas,imol,jmol,npx,npy,npz,nptbas
      integer j,ntot,npairs,j1,j2,j3,j4,l,first,last,maxatoms
      integer multiplicity
      integer, dimension (:), allocatable :: val_label
      integer*2 ind,jnd
      character*16 atdum
      character*1 molecule(40), yes, blank,atom(16)
      character*4 end_data
      character*40 dumname, filebas,filrbas,filegam,filesyn
      integer endx,endy,endz,endij,ij,il,jk,ik,jl,munu,nbas2
      integer i,mu,nu,ix,iy,iz,kmol,lmol,kl,ki,lj
      integer n_elf, n_aim, atnum(5)
      logical jlimit,if_file,elf_pop,rho_pop,elf_and_rho,if_calc
      common/grid/ dx,dy,dz
      common /files/ molecule, filebas
      common /erftab/ erf0,erf1
      equivalence (molecule(1),dumname),(atom(1),atdum)
      data blank /' '/
      open(unit=45,file='elf09_bas.gjf',status='old',access='append')
*     read wave function information
      call read_wfn
      multiplicity=n_alpha-n_beta+1
      nmola=nmol
      maxatoms=natoms+10
      maxbasin=400
      allocate(vol_at(maxatoms),yij(maxatoms,maxatoms))
      allocate(cc(nmol*(nmol+1)/2),rho_at(maxatoms))
      allocate(taa(maxatoms),tbb(maxatoms),tab(maxatoms),tbas(maxatoms))
      allocate(t_i(maxatoms),ta(maxatoms),tb(maxbasin))
      allocate(t_ij(maxatoms,maxatoms),aim_ind(maxatoms))
      allocate(indice_ij(nmol*(nmol+1)/2))
      allocate(aaij(maxbasin,maxbasin),abij(maxbasin,maxbasin))
      allocate(bbij(maxbasin,maxbasin))
      allocate(indij(nmol*(nmol+1)/2),indi(nmol*(nmol+1)/2))
      allocate(smunu(ngto*(ngto+1)/2))
      allocate(tij(maxatoms,nmol*(nmol+1)/2))
      allocate(name_atom(maxatoms))
      do i=1,nmol 
         indij(i)=(i-1)*(nmol+nmol-i)/2
         do j=imol,nmol
            k=indij(i)+i
            indi(k)=i
         enddo
      enddo
      allocate(induhfa(nmol*(nmol+1)/2))
      allocate(induhfb(nmol*(nmol+1)/2))
      do i=1,n_alpha
         do j=i,n_alpha
           ij=(i-1)*(nmol+nmol-i)/2+j
           induhfa(ij)=indij(i)+j
           induhfb(ij)=indij(i+n_alpha)+j+n_alpha
         enddo
      enddo
*      write(*,*) ' threshold for integration tol ==> 10-tol '
*      read(*,*) tol
       tol=7
      write(*,*) 'Number of ELF and AIM basins considered'
      read(*,*)  n_elf, n_aim
*        n_elf = 0
*        n_aim = 0
      if(n_elf.ne.0) then
         allocate(elf_ind(n_elf))
         write(*,*) 'ELF basins considered '
         read(*,*) (elf_ind(i),i=1,n_elf)
      endif
      if(n_aim.ne.0) then
*         allocate(aim_ind(n_aim))
         write(*,*) 'AIM basins considered '
         read(*,*) (aim_ind(i),i=1,n_aim)
      endif
*
*     create input file name
*
      i=1
      do while (molecule(i).eq.blank)
         do j=i,39
            molecule(j)=molecule(j+1)
         enddo
         i=i+1
      enddo
      do while (molecule(i).ne.blank)
         i=i+1
      enddo
      molecule(i)='_'
      molecule(i+1)='e'
      molecule(i+2)='b'
      molecule(i+3)='a'
      molecule(i+4)='s'
      molecule(i+5)='.'
      molecule(i+6)='s'
      molecule(i+7)='b'
      molecule(i+8)='f'
      do j=i+9,40
         molecule(j)=blank
      enddo
      jnd=i
      filebas=dumname
      molecule(i+1)='r'
      filrbas=dumname
      molecule(i+1)='g'
      molecule(i+2)='a'
      molecule(i+3)='m'
      molecule(i+4)='.'
      molecule(i+5)='s'
      molecule(i+6)='b'
      molecule(i+7)='f'
      do j=i+8,40
         molecule(j)=blank
      enddo
      filegam=dumname
      molecule(i+1)='e'
      molecule(i+2)='s'
      molecule(i+3)='y'
      molecule(i+4)='n'
      molecule(i+5)='.'
      molecule(i+6)='s'
      molecule(i+7)='b'
      molecule(i+8)='f'
      do j=i+9,40
         molecule(j)=blank
      enddo
      filesyn=dumname
*
*     input localization basin indexes from sbf file 
*
      elf_pop=.true.
      elf_and_rho=.true.
      write(*,'(2x,a40)') filebas
      inquire(file=filebas,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist: ",a40)') filebas
        elf_pop=.false.
        elf_and_rho=.false.
      endif
      if(elf_pop) then
         open(unit=1,file=filebas,status='old',form='unformatted')
         read(1) npx,npy,npz
         np(1)=npx
         np(2)=npy
         np(3)=npz
         read(1) xmin,xmax,ymin,ymax,zmin,zmax
         npts=npx*npy*npz
         nptbas=npts
         allocate(attract_code(npts))
         allocate(ind_at(npts))
         read(1) k,nbas
         allocate(val_label(nbas))
         allocate(sij(nbas,nmol*(nmol+1)/2))
         allocate(name_basin(nbas),volume(nbas),x(nbas),y(nbas),z(nbas))
         allocate(temp_label(nbas),bij(nbas,nbas))
         allocate(stij(nbas,nbas,nmol),rho_bas(nbas),sbas(nbas))
         allocate(p_i(nbas),p_ij(nbas,nbas),stbas(nbas,maxatoms))
         allocate(at_contrib(nbas,maxatoms))
         allocate(gx(npx),gy(npy),gz(npz))
         allocate(indice_x(npx),indice_y(npy),indice_z(npz))
         allocate(popa(nbas),popb(nbas),popaa(nbas),popbb(nbas))
         allocate(stdev(nbas),popab(nbas))
         allocate(cova(nbas,nbas+1),covb(nbas,nbas),covab(nbas,nbas+1))
         do i=1,k
            read(1) nat(i),end_data,weight,weight,weight
         enddo
         nbas2=0
         do i=1,nbas
            read(1) val_label(i),name_basin(i),volume(i),x(i),y(i),z(i)
         enddo
         do i=1,nbas
            if_calc=.false.
            if(n_elf.gt.0) then
               do j=1,n_elf
                  if(elf_ind(j).eq.i) if_calc=.true.
               enddo
            else
               if_calc=.true.
            endif
            if(volume(i).le.1.0d-06) if_calc=.false.
            if(if_calc) then
               nbas2=nbas2+1
               temp_label(i)=nbas2
               name_basin(nbas2)=name_basin(i)
               volume(nbas2)=volume(i)
               x(nbas2)=x(i)
               y(nbas2)=y(i)
               z(nbas2)=z(i)
            else
               temp_label(i)=0
            endif
         enddo
         read(1)(attract_code(i),i=1,npts)
         close(unit=1)
         do i=1,npts
            k=attract_code(i)
            if(k.gt.0) then
               attract_code(i)=temp_label(k)
            endif
         enddo
         nbas=nbas2
         write(*,'(3x,'' number of localization basins  '',i3)') nbas
      else
         do i=1,npts
            attract_code(i)=0
         enddo
      endif
      rho_pop=.true.
      write(*,'(2x,a40)') filrbas
      inquire(file=filrbas,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist: ",a40)') filrbas
        rho_pop=.false.
        elf_and_rho=.false.
      endif
      if(rho_pop) then
         open(unit=1,file=filrbas,status='old',form='unformatted')
         read(1) npx,npy,npz
         read(1) xmin,xmax,ymin,ymax,zmin,zmax
         npts=npx*npy*npz
         if(npts.ne.nptbas) then
            write(*,*) 'elf and rho bas files have different size '
            stop
         endif
         read(1) k,natoms
c        allocate(atom_name(natoms))
         do i=1,k
            read(1) nat(i),atom_name(i),weight,weight,weight
         enddo
         nbas2=0
         do i=1,natoms
            read(1) k,name_atom(i),vol_at(i),xat(i),yat(i),zat(i)
            if_calc=.false.
            if(n_aim.gt.0) then
               do j=1,n_aim
                  if(aim_ind(j).eq.i) if_calc=.true.
               enddo
            else
               if_calc=.true.
            endif
            if(vol_at(i).le.1.0d-06) if_calc=.false.
            if(if_calc) then
               nbas2=nbas2+1
               temp_label(i)=nbas2
               name_atom(nbas2)=name_atom(i)
               vol_at(nbas2)=vol_at(i)
               x(nbas2)=x(i)
               y(nbas2)=y(i)
               z(nbas2)=z(i)
            else
               temp_label(i)=0
            endif
         enddo
         read(1)(ind_at(i),i=1,npts)
         j=0
         do i=1,npts
            k=ind_at(i)
            if(k.gt.0) then
               ind_at(i)=temp_label(k)
               j=j+1
            endif
         enddo
         close(unit=1)
         natoms=nbas2
         write(*,'(3x,'' number of atomic basins  '',i3)') natoms
      else
         if(elf_pop) then
            do i=1,npts
               ind_at(i)=0
            enddo
         else
            write(*,*) ' no bas.sbf datasets stop '
            stop
         endif
      endif
      do i=1,natoms
         atdum=name_atom(i)
         do j=1,4
            atom(j)=blank
         enddo
         do j=6,9
            if(atom(j).ne.')') then
               atom(j-5)=atom(j)
            endif
         enddo   
         name_atom(i)=atom(1)//atom(2)//atom(3)//atom(4)
      enddo
*     read density matrices coefficients
      inquire(file=filegam,exist=if_file)
      if(if_file) then
      write(*,'(2x,"correlated wave function coeff ",a40)') filegam
      open(unit=14,file=filegam,status='unknown',form='unformatted')
      read(14) na
      allocate(ia(na),pa(na))
      read(14)(pa(i),i=1,na)
      read(14)(ia(i),i=1,na)
      if(uhf) then
         do i=1,na
            j=ia(i)
            ia(i)=induhfa(j)
         enddo
      endif
      read(14) nb
      allocate(ib(nb),pb(nb))
      nmola=nmol
      if(nb.ne.0) then
         nmola=nmol/2
         read(14)(pb(i),i=1,nb)
         read(14)(ib(i),i=1,nb)
      endif
      if(uhf) then
         do i=1,nb
            j=ib(i)
            ib(i)=induhfb(j)
         enddo
      endif
*     input 2nd density matrix coefficients
      read(14) naa
      allocate(iaa(naa),jaa(naa),kaa(naa),laa(naa),paa(naa))
      read(14)(paa(i),i=1,naa)
      read(14)(iaa(i),i=1,naa)
      read(14)(jaa(i),i=1,naa)
      read(14)(kaa(i),i=1,naa)
      read(14)(laa(i),i=1,naa)
      read(14) nbb
      if(nbb.ne.0) then
         allocate(ibb(nbb),jbb(nbb),kbb(nbb),lbb(nbb),pbb(nbb))
         read(14)(pbb(i),i=1,nbb)
         read(14)(ibb(i),i=1,nbb)
         read(14)(jbb(i),i=1,nbb)
         read(14)(kbb(i),i=1,nbb)
         read(14)(lbb(i),i=1,nbb)
         if(uhf) then
            do i=1,nbb
               ibb(i)=ibb(i)+n_alpha
               jbb(i)=jbb(i)+n_alpha
               kbb(i)=kbb(i)+n_alpha
               lbb(i)=lbb(i)+n_alpha
            enddo
         endif
      endif
      read(14) nab
      allocate(iab(2*nab),jab(2*nab),kab(2*nab),lab(2*nab),pab(2*nab))
      read(14)(pab(i),i=1,nab)
      read(14)(iab(i),i=1,nab)
      read(14)(jab(i),i=1,nab)
      read(14)(kab(i),i=1,nab)
      read(14)(lab(i),i=1,nab)
      ifno=.false.
      close(unit=14)
      else if(open_shell) then
         allocate(iaa(nmol*(nmol+1)/2),jaa(nmol*(nmol+1)/2))
         allocate(kaa(nmol*(nmol+1)/2),laa(nmol*(nmol+1)/2))
         allocate(ibb(nmol*(nmol+1)/2),jbb(nmol*(nmol+1)/2))
         allocate(kbb(nmol*(nmol+1)/2),lbb(nmol*(nmol+1)/2))
         allocate(iab(nmol*(nmol+1)),jab(nmol*(nmol+1)))
         allocate(kab(nmol*(nmol+1)),lab(nmol*(nmol+1)))
         allocate(paa(nmol*(nmol+1)/2),pbb(nmol*(nmol+1)/2))
         allocate(pab(nmol*(nmol+1)))
         if(uhf) then
            na=0
            nb=0
            naa=0
            nbb=0
            nab=0
            nmola=n_alpha
            allocate(ia(n_alpha),pa(n_alpha),ib(n_beta),pb(n_beta))
            do imol=1,n_alpha
               if(occ(1,imol).gt.0.0d0) then
                  na=na+1
                  pa(na)=1.0d0
                  if(ifno) then
                    pa(na)=occ(1,imol)
                  endif
                  ia(na)=indij(imol)+imol
                  do jmol=imol,n_alpha
                     if(occ(1,jmol).gt.0.0d0) then
                        naa=naa+1
                        paa(naa)=1.0d0
                        iaa(naa)=imol
                        jaa(naa)=jmol
                        kaa(naa)=imol
                        laa(naa)=jmol
                        if(ifno) then
                           tota=occ(1,jmol)*occ(1,imol)
                           paa(naa)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
            do imol=n_alpha+1,n_alpha+n_beta
               if(occ(2,imol).gt.0.0d0) then
                  nb=nb+1
                  pb(nb)=1.0d0
                  if(ifno) then
                     pb(nb)=occ(2,imol)
                  endif
                  ib(nb)=indij(imol)+imol
                  do jmol=imol,n_alpha+n_beta
                     if(occ(2,jmol).gt.0.0d0) then
                        nbb=nbb+1
                        pbb(nbb)=1.0d0
                        ibb(nbb)=imol
                        jbb(nbb)=jmol
                        kbb(nbb)=imol
                        lbb(nbb)=jmol
                        if(ifno) then
                           tota=occ(2,jmol)*occ(2,imol)
                           pbb(nbb)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
            do imol=1,n_alpha
               do jmol=n_alpha+1,n_alpha+n_beta
                  tota=occ(1,imol)*occ(2,jmol)
                  if(tota.gt.1.0d-6) then
                     nab=nab+1
                     iab(nab)=imol
                     jab(nab)=jmol
                     kab(nab)=imol
                     lab(nab)=jmol
                     pab(nab)=tota
                  endif
               enddo
            enddo
****     end UHF 
         else
            na=0
            nb=0
            naa=0
            nbb=0
            nmola=nmol
            allocate(ia(nmol),pa(nmol),ib(nmol),pb(nmol))
            do imol=1,nmol
               if(occ(1,imol).gt.0.0d0) then
                  na=na+1
                  pa(na)=1.0d0
                  if(ifno) then
                     pa(na)=occ(1,imol)
                  endif
                  ia(na)=indij(imol)+imol
                  do jmol=imol,nmol
                     if(occ(1,jmol).gt.0.0d0) then
                        naa=naa+1
                        paa(naa)=1.0d0
                        iaa(naa)=imol
                        jaa(naa)=jmol
                        kaa(naa)=imol
                        laa(naa)=jmol
                        if(ifno) then
                           tota=occ(1,jmol)*occ(1,imol)
                           paa(naa)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
            do imol=1,nmol
               if(occ(2,imol).gt.0.0d0) then
                  nb=nb+1
                  pb(nb)=1.0d0
                  if(ifno) then
                     pb(nb)=occ(2,imol)
                  endif
                  ib(nb)=indij(imol)+imol
                  do jmol=imol,nmol
                     if(occ(2,jmol).gt.0.0d0) then
                        nbb=nbb+1
                        pbb(nbb)=1.0d0
                        ibb(nbb)=imol
                        jbb(nbb)=jmol
                        kbb(nbb)=imol
                        lbb(nbb)=jmol
                        if(ifno) then
                           tota=occ(2,jmol)*occ(2,imol)
                           pbb(nbb)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
            do imol=1,nmol
               do jmol=imol,nmol
                  tota=occ(1,imol)*occ(2,jmol)+occ(2,imol)*occ(1,jmol)
                  if(tota.gt.1.0d-6) then
                     nab=nab+1
                     iab(nab)=imol
                     jab(nab)=jmol
                     kab(nab)=imol
                     lab(nab)=jmol
                     pab(nab)=0.5d0*tota
                  endif
               enddo
            enddo
         endif
***      end open shell
      else
        na=0
        nb=0
        naa=0
        nbb=0
        nab=0
        allocate(ia(nmol),pa(nmol),ib(nmol),pb(nmol))
        allocate(iaa(nmol*(nmol+1)/2),jaa(nmol*(nmol+1)/2))
        allocate(kaa(nmol*(nmol+1)/2),laa(nmol*(nmol+1)/2))
        allocate(iab(nmol*(nmol+1)/2),jab(nmol*(nmol+1)/2))
        allocate(kab(nmol*(nmol+1)/2),lab(nmol*(nmol+1)/2))
        allocate(paa(nmol*(nmol+1)/2),pab(nmol*(nmol+1)/2))
        do imol=1,nmol
           if(occ(1,imol).gt.0.0d0) then
              na=na+1
              pa(na)=1.0d0
              ia(na)=indij(imol)+imol
              if(ifno) then
                 pa(na)=occ(1,imol)
              endif
              do jmol=imol,nmol
                 if(occ(1,jmol).gt.0.0d0) then
                    tota=occ(1,jmol)*occ(2,imol)
                    naa=naa+1
                    paa(naa)=1.0d0
                    iaa(naa)=imol
                    jaa(naa)=jmol
                    kaa(naa)=imol
                    laa(naa)=jmol
                    if(ifno) then
                       paa(naa)=tota
                    endif
                    nab=nab+1
                    pab(nab)=1.0d0
                    iab(nab)=imol
                    jab(nab)=jmol
                    kab(nab)=imol
                    lab(nab)=jmol
                    if(ifno) then
                       paa(naa)=tota
                       pab(nab)=tota
                    endif
                 endif
              enddo
           endif
        enddo
      endif
      if(open_shell) then
         if(.not.uhf) then
           if(ifno) then
              ifno=.false.
           endif
         endif
      endif
      call xtime(t)
      spent=t(1)+t(2)
      rho_tol=0.0d0
      pi=dacos(-1.0d0)
      const=2.0d0/dsqrt(pi)
      do j=0,10000
         arg=5.0d-4*dfloat(j)
         erf0(j)=derf(arg)
         erf1(j)=const*dexp(-arg*arg)
      enddo
      total_charge=0.0d0
      do i=1,n_alpha
         total_charge=total_charge+occ(1,i)
      enddo
      do i=1,n_beta
         total_charge=total_charge+occ(2,i)
      enddo
      call box_dens(rho_tol)
      munu=0
      do mu=1,ngto
      do nu=mu,ngto
         munu=munu+1
      enddo
      enddo
      do imol=1,nmol
         do jmol=imol,nmol
            lmol=indij(imol)+jmol
            munu=0
            sij(1,lmol)=0.0d0
            do mu=1,ngto
               do nu=mu,ngto
                  munu=munu+1
                  if(mu.eq.nu) then
                     arg=c(imol,mu)*c(jmol,nu)*smunu(munu)
                  else
                     arg=(c(imol,mu)*c(jmol,nu)+
     .                    c(jmol,mu)*c(imol,nu))*smunu(munu)
                  endif
                  sij(1,lmol)=sij(1,lmol)+arg
               enddo
            enddo
            if(imol.ne.jmol) then
               sij(1,lmol)=sij(1,lmol)+sij(1,lmol)
            endif
         enddo
      enddo
*     alpha contributions
      totala=0.0d0
      do l=1,na
         lmol=ia(l)
         totala=totala+pa(l)*sij(1,lmol)
      enddo
      totalaa=0.0d0
      do l=1,naa
         imol=iaa(l)
         jmol=jaa(l)
         kmol=kaa(l)
         lmol=laa(l)
         amult=2.0d0
         ik=indij(imol)+kmol
         if(imol.gt.kmol) then
            ik=indij(kmol)+imol
         endif
         jl=indij(jmol)+lmol
         if(jmol.gt.lmol) then
            jl=indij(lmol)+jmol
         endif
         il=indij(imol)+lmol
         if(imol.gt.lmol) then
            il=indij(lmol)+imol
         endif
         jk=indij(jmol)+kmol
         if(jmol.gt.kmol) then
            jk=indij(kmol)+jmol
         endif
         if(ik.eq.jl) then
            amult=1.0d0
         endif
         if(ifno) then
            totalaa=totalaa+amult*(paa(l)*sij(1,ik)*sij(1,jl)-
     .                   dsqrt(paa(l))*sij(1,il)*sij(1,jk))
         else
            totalaa=totalaa+amult*paa(l)*(sij(1,ik)*sij(1,jl)-
     .                                 sij(1,il)*sij(1,jk))
         endif
      enddo
*     sum for beta spin
      if(nb.ne.0) then
      totalb=0.0d0
      do l=1,nb
         lmol=ib(l)
         totalb=totalb+pb(l)*sij(1,lmol)
      enddo
      totalbb=0.0d0
      do l=1,nbb
         imol=ibb(l)
         jmol=jbb(l)
         kmol=kbb(l)
         lmol=lbb(l)
         amult=2.0d0
         ik=indij(imol)+kmol
         if(imol.gt.kmol) then
            ik=indij(kmol)+imol
         endif
         jl=indij(jmol)+lmol
         if(jmol.gt.lmol) then
            jl=indij(lmol)+jmol
         endif
         il=indij(imol)+lmol
         if(imol.gt.lmol) then
            il=indij(lmol)+imol
         endif
         jk=indij(jmol)+kmol
         if(jmol.gt.kmol) then
            jk=indij(kmol)+jmol
         endif
         if(ik.eq.jl) then
            amult=1.0d0
         endif
         if(ifno) then
            totalbb=totalbb+amult*(pbb(l)*sij(1,ik)*sij(1,jl)-
     .                    dsqrt(pbb(l))*sij(1,il)*sij(1,jk))
         else
            totalbb=totalbb+amult*pbb(l)*(sij(1,ik)*sij(1,jl)-
     .                                 sij(1,il)*sij(1,jk))
         endif
      enddo
      else
         totalb=totala
         totalbb=totalaa
      endif
*     sum for opposite spin
      do imol=1,nmol-1
         do jmol=imol+1,nmol
            lmol=indij(imol)+jmol
         enddo
      enddo
      totalab=0.0d0
      do l=1,nab
         imol=iab(l)
         jmol=jab(l)
         kmol=kab(l)
         lmol=lab(l)
         amult=2.0d0
         ik=indij(imol)+kmol
         if(imol.gt.kmol) then
            ik=indij(kmol)+imol
         endif
         jl=indij(jmol)+lmol
         if(jmol.gt.lmol) then
            jl=indij(lmol)+jmol
         endif
         if(uhf.and.ifci) then
            jl=indij(jmol+n_alpha)+lmol+n_alpha
         endif
         if(ik.eq.jl) then
            amult=1.0d0
         endif
         totalab=totalab+amult*pab(l)*sij(1,ik)*sij(1,jl)
      enddo
      if(.not.uhf) then
      totalab=2.0d0*totalab
      endif
      units=1.0d0
      total=totala+totalb
      arg=totalaa+totalbb+totalab+total*(1.0d0-total)
      write(*,'(3x,'' total integrated density '',f12.6)') total
      write(*,'(3x,'' total variance '',f12.6)') arg
      if((total_charge-total).gt.0.2d0) then
         write(*,*) ' ** WARNING increase the box size '
      endif
      if(open_shell) then
         totals=0.5d0*(totala-totalb)
         write(*,'(3x,'' total integrated spin density '',f12.6)') 
     .                                                      totals
      endif
      ntot=dint(total)
      npairs=ntot*(ntot-1)
      tol=1.0d0/1.0d1**tol
      dx=(xmax-xmin)/(npx-1)
      dy=(ymax-ymin)/(npy-1)
      dz=(zmax-zmin)/(npz-1)
      if(elf_pop) then
         do i=1,nbas
            do j=1,(nmol*(nmol+1))/2
               sij(i,j)=0.0d0
            enddo
            do j=1,nmol
               if(elf_and_rho) then
                  do k=1,natoms
                     stij(i,k,j)=0.0d0
                  enddo
               endif
            enddo
            rho_bas(i)=0.0d0
         enddo
      endif
      if(rho_pop) then
         do i=1,natoms
            do j=1,(nmol*(nmol+1))/2
               tij(i,j)=0.0d0
            enddo
            rho_at(i)=0.0d0
         enddo
      endif
      munu=0
      do mu=1,ngto
         do nu=mu,ngto
            munu=munu+1
            if(elf_pop) then
               do i=1,nbas
                  sbas(i)=0.0d0
                  if(elf_and_rho) then
                     do j=1,natoms
                        stbas(i,j)=0.0d0
                     enddo
                  endif
               enddo
            endif
            if(rho_pop) then
               do i=1,natoms
                  tbas(i)=0.0d0
               enddo
            endif
            call munu_contrib(mu,nu)
            endx=0
            do i=1,npx
               if(dabs(gx(i)).ge.tol) then
                 endx=endx+1
                 indice_x(endx)=i
               endif
            enddo
            endy=0
            do i=1,npy
               if(dabs(gy(i)).ge.tol) then
                 endy=endy+1
                 indice_y(endy)=i
               endif
            enddo
            endz=0
            do i=1,npz
               if(dabs(gz(i)).ge.tol) then
                 endz=endz+1
                 indice_z(endz)=i
               endif
            enddo
            do i=1,(nmol*(nmol+1))/2
               cc(ij)=0.0d0
            enddo
            if(mu.eq.nu) then
              do i=1,nmol
                 do j=i,nmol
                    ij=indij(i)+j
                    cc(ij)=c(i,mu)*c(j,mu)
                 enddo
              enddo
            else
              do i=1,nmol
                 do j=i,nmol
                    ij=indij(i)+j
                    cc(ij)=c(i,mu)*c(j,nu)+c(i,nu)*c(j,mu)
                 enddo
              enddo
            endif
            do i=1,endx
               ix=indice_x(i)
               do j=1,endy
                  iy=indice_y(j)
                  do k=1,endz
                     iz=indice_z(k)
                     npt=iz+npz*((iy-1)+npy*(ix-1))
                     ind=attract_code(npt)
                     jnd=ind_at(npt)
                     rho_box=gx(ix)*gy(iy)*gz(iz)
                     if(ind.gt.0) then
                        sbas(ind)=sbas(ind)+rho_box
                        if(jnd.gt.0) then
                           stbas(ind,jnd)=stbas(ind,jnd)+rho_box
                        endif
                     endif
                     if(jnd.gt.0) then
                        tbas(jnd)=tbas(jnd)+rho_box
                     endif
                  enddo
               enddo
            enddo
            do ind=1,nbas
               do imol=1,nmol
                  ij=indij(imol)+imol
                  do jnd=1,natoms
                     stij(ind,jnd,imol)=stij(ind,jnd,imol)+
     .                                  cc(ij)*stbas(ind,jnd)
                  enddo
                  do jmol=imol,nmol
                     ij=indij(imol)+jmol
                     sij(ind,ij)=sij(ind,ij)+cc(ij)*sbas(ind)
                  enddo
               enddo
            enddo
            do ind=1,natoms
               do imol=1,nmol
                  do jmol=imol,nmol
                     ij=indij(imol)+jmol
                     tij(ind,ij)=tij(ind,ij)+cc(ij)*tbas(ind)
                  enddo
               enddo
            enddo
         enddo
      enddo
      if(elf_pop) then
      total_charge=total
      do i=1,nbas
         totala=0.0d0
         totalb=0.0d0
         totalaa=0.0d0
         totalbb=0.0d0
         totalab=0.0d0
         do j=1,natoms
            at_contrib(i,j)=0.0d0
         enddo
         do l=1,na
            lmol=ia(l)
            totala=totala+pa(l)*sij(i,lmol)
            k=indi(lmol)
            do j=1,natoms
               at_contrib(i,j)=at_contrib(i,j)+pa(l)*stij(i,j,k)
            enddo
         enddo
         do l=1,naa
            imol=iaa(l)
            jmol=jaa(l)
            kmol=kaa(l)
            lmol=laa(l)
            amult=2.0d0
            ik=indij(imol)+kmol
            if(imol.gt.kmol) then
               ik=indij(kmol)+imol
            endif
            jl=indij(jmol)+lmol
            if(jmol.gt.lmol) then
               jl=indij(lmol)+jmol
            endif
            il=indij(imol)+lmol
            if(imol.gt.lmol) then
               il=indij(lmol)+imol
            endif
            jk=indij(jmol)+kmol
            if(jmol.gt.kmol) then
               jk=indij(kmol)+jmol
            endif
            if((i.eq.j).and.(k.eq.l)) then
               amult=1.0d0
            endif
            if(ifno) then
               totalaa=totalaa+amult*(paa(l)*sij(i,ik)*sij(i,jl)-
     .                       dsqrt(paa(l))*sij(i,il)*sij(i,jk))
            else
               totalaa=totalaa+amult*paa(l)*(sij(i,ik)*sij(i,jl)-
     .                                 sij(i,il)*sij(i,jk))
            endif
         enddo
         if(nb.ne.0) then
            do l=1,nb
               lmol=ib(l)
               totalb=totalb+pb(l)*sij(i,lmol)
               k=indi(lmol)
               do j=1,natoms
                  at_contrib(i,j)=at_contrib(i,j)+pb(l)*stij(i,j,k)
               enddo
            enddo
            totalbb=0.0d0
            do l=1,nbb
               imol=ibb(l)
               jmol=jbb(l)
               kmol=kbb(l)
               lmol=lbb(l)
               amult=2.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=1.0d0
               endif
               if(ifno) then
                  totalbb=totalbb+amult*(pbb(l)*sij(i,ik)*sij(i,jl)-
     .                         dsqrt(pbb(l))*sij(i,il)*sij(i,jk))
               else
                  totalbb=totalbb+amult*pbb(l)*(sij(i,ik)*sij(i,jl)-
     .                                 sij(i,il)*sij(i,jk))
               endif
            enddo
         else
            totalb=totala
            totalbb=totalaa
            do j=1,natoms
               at_contrib(i,j)=2.0d0*at_contrib(i,j)
            enddo
         endif
*     sum for opposite spin
         do l=1,nab
            imol=iab(l)
            jmol=jab(l)
            kmol=kab(l)
            lmol=lab(l)
            amult=2.0d0
            ik=indij(imol)+kmol
            if(imol.gt.kmol) then
               ik=indij(kmol)+imol
            endif
            jl=indij(jmol)+lmol
            if(jmol.gt.lmol) then
               jl=indij(lmol)+jmol
            endif
            if(uhf.and.ifci) then
               jl=indij(jmol+n_alpha)+lmol+n_alpha
            endif
            if(ik.eq.jl) then
               amult=1.0d0
            endif
            totalab=totalab+amult*pab(l)*sij(i,ik)*sij(i,jl)
         enddo
         total=totala+totalb
         rho_bas(i)=total
         popa(i)=totala
         popb(i)=totalb
         if(uhf) then
            popab(i)=totalab
         else
            popab(i)=2.0d0*totalab
         endif
         popaa(i)=totalaa
         popbb(i)=totalbb
         arg=popaa(i)+popbb(i)+popab(i)+rho_bas(i)*(1.0d0-rho_bas(i))
         bij(i,i)=arg
         stdev(i)=dsqrt(arg)
         if(rho_bas(i).eq.0.0d0) then
           rho_bas(i)=1.0d-3
         endif
      enddo
***** rename ELF basins
      if(elf_and_rho) then
         do i=1,nbas
            if(val_label(i).ne.5) then
               val_label(i)=0
               k=0
               do j=1,natoms
                  weight=at_contrib(i,j)/rho_bas(i)
                  if(weight.ge.0.05d0) then
                     val_label(i)=val_label(i)+1
                     atnum(val_label(i))=j
                     if(nat(j).eq.1) then
                        k=2
                     endif
                  endif
                  if(val_label(i).eq.1) then
                     dumname='V('//atom_name(atnum(1))//')'
                  else if (val_label(i).eq.2) then
                     dumname='V('//atom_name(atnum(1))//','//
     .                             atom_name(atnum(2))//')'
                  else if (val_label(i).ge.3) then
                     dumname='V('//atom_name(atnum(1))//','//
     .                             atom_name(atnum(2))//','//
     .                             atom_name(atnum(3))//')'
                  endif
                  call short_name(dumname,name_basin(i))
               enddo
               if(val_label(i).gt.1) then
                  val_label(i)=3
                  if(k.eq.2) val_label(i)=2
               endif
            endif
         enddo
         do i=1,npts
            k=attract_code(i)
            if(k.ne.0) then
               attract_code(i)=val_label(k)
            else
               attract_code(i)=0
            endif
         enddo
         open(unit=2,file=filesyn,status='unknown',form='unformatted')
         write(2) np
         write(2) xmin,xmax,ymin,ymax,zmin,zmax
         write(2) (attract_code(i),i=1,npts)
         close(unit=2)
      endif
      total=0.0d0
      write(45,*) 'Basin populations'
      write(45,'(//,t3,''basin'',t23,''vol.'',t31,''pop.'',t39,"pab",
     .          t47,"paa",t54,"pbb",t62,''sigma2'',
     .          t70,''std. dev.'',//)')
      write(*,'(//,t3,''basin'',t23,''vol.'',t31,''pop.'',t39,"pab",
     .          t47,"paa",t54,"pbb",t62,''sigma2'',
     .          t70,''std. dev.'',//)')
      do i=1,nbas
         write(*,'(i3,1x,a16,f6.2,6(2x,f6.2))') i,name_basin(i),
     .                 volume(i),rho_bas(i),popab(i),popaa(i),popbb(i),
     .                 bij(i,i),stdev(i)
         total=total+rho_bas(i)
         write(45,'(i3,1x,a16,f6.2,6(2x,f6.2))') i,name_basin(i),
     .                 volume(i),rho_bas(i),popab(i),popaa(i),popbb(i),
     .                 bij(i,i),stdev(i)
      enddo
      write(*,'(//,'' sum of populations '',f12.6,/)') total
      if(dabs(total_charge-total).ge.0.2d0) then
         write(*,*) ' ** WARNING increase integration threshold '
         stop
      endif
      write(*,'(//,t3,''basin'',t23,''pa'',t31,''paa.'',t39,"sigma2",
     .          t47,"pb",t54,"pbb",t62,''sigma2'',//)')
      do i=1,nbas
         arg=popaa(i)+popa(i)*(1.0d0-popa(i))
         total=popbb(i)+popb(i)*(1.0d0-popb(i))
         write(*,'(i3,1x,a16,6(f6.2,2x))') i,name_basin(i),popa(i),
     .                                     popaa(i),
     .                                     arg,popb(i),popbb(i),total
      enddo
      if(.not.ifci) then
      write(*,'(//,'' orbital contributions '',/)')
      j1=-9
      j2=0
      jlimit=.true.
      imol=nmol
      if(fukui) imol=nmol-1
      allocate(pop_orb(nmol))
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nmol) then
            j2=nmol
            jlimit=.false.
         endif
         write(*,'(//,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            if(j2.le.imol) then
               do j=j1,j2
                  ij=nmol*(j-1)+(3*j-j*j)/2
                  pop_orb(j)=(occ(1,j)+occ(2,j))*sij(i,ij)
               enddo
            else
               ij=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
               pop_orb(nmol)=sij(i,ij)
            endif
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                                (pop_orb(j),j=j1,j2)
        enddo
      enddo
      if(fukui) then
         if(open_shell) then
         write(*,'(//,'' condensed Fukui indexes '',/)')
         j1=nmol*(imol-1)+(3*imol-imol*imol)/2
         j2=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
         j3=nmol*(n_beta-1)+(3*n_beta-n_beta*n_beta)/2
         j4=nmol*n_beta+(n_beta-n_beta**2+2)/2
         if(multiplicity.gt.2) then
            write(*,'(t22,''f-'',t28,''f+'',t34,"fs+",t39,"fs-",//)')
            do i=1,nbas
            write(*,'(1x,a16,4f6.2)')name_basin(i),sij(i,j1),sij(i,j4),
     .                                0.5d0*(sij(i,j2)+sij(i,j3)),
     .                                0.5d0*(sij(i,j1)+sij(i,j4))
            enddo
            else
            write(*,'(t22,''f-'',t28,''f+'',t34,"fs+",t39,//)')
            do i=1,nbas
            write(*,'(1x,a16,3f6.2)')name_basin(i),sij(i,j1),sij(i,j4),
     .                                0.5d0*(sij(i,j2)+sij(i,j3))
            enddo
         endif
         else
         write(*,'(//,'' condensed Fukui indexes '',/)')
         j1=nmol*(imol-1)+(3*imol-imol*imol)/2
         j2=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
         write(*,'(t22,''f-'',t28,''f+'',t34,''f0'',t39,''Df'',//)')
         do i=1,nbas
            write(*,'(1x,a16,4f6.2)')name_basin(i),sij(i,j1),sij(i,j2),
     .                                0.5d0*(sij(i,j1)+sij(i,j2)),
     .                                sij(i,j1)-sij(i,j2)
         enddo
         endif
      endif
      endif
      if(elf_and_rho) then
      write(*,'(//,'' atomic contributions '',/)')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(//,18x,10(1x,a4,1x))')(name_atom(j),j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                                (at_contrib(i,j),j=j1,j2)
        enddo
      enddo
      endif
      if(open_shell) then
         write(*,'(//,t1,''basin'',t20,''integrated spin density'',//)'
     .            )
         totals=0.0d0
         do i=1,nbas
            total=popa(i)-popb(i)
            total=0.5d0*total
            write(*,'(1x,a16,f12.6)') name_basin(i),total
            totals=totals+total
         enddo
         write(*,'(//,'' sum of spin densities '',f12.6)') totals
      endif
      do i=1,nbas
         do j=1,nbas
            totalaa=0.0d0
            do l=1,naa
               imol=iaa(l)
               jmol=jaa(l)
               kmol=kaa(l)
               lmol=laa(l)
               amult=1.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=0.5d0
               endif
               if(ifno) then
                  totalaa=totalaa+amult*(paa(l)*(sij(i,ik)*sij(j,jl)+
     .                  sij(j,ik)*sij(i,jl))-dsqrt(paa(l))*(sij(i,il)*
     .                  sij(j,jk)+sij(j,il)*sij(i,jk)))
               else
                  totalaa=totalaa+amult*paa(l)*(sij(i,ik)*sij(j,jl)+
     .                   sij(j,ik)*sij(i,jl)-sij(i,il)*sij(j,jk)-
     .                   sij(j,il)*sij(i,jk))
               endif
            enddo
            cova(i,j)=totalaa-popa(i)*popa(j)
            aaij(i,j)=totalaa
         enddo
         cova(i,i)=cova(i,i)+popa(i)
         stdev(i)=dsqrt(cova(i,i))
      enddo
      write(*,'(//,'' alpha spin covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                                      (cova(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                           (cova(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      if(open_shell) then
      do i=1,nbas
         do j=1,nbas
            totalbb=0.0d0
            do l=1,nbb
               imol=ibb(l)
               jmol=jbb(l)
               kmol=kbb(l)
               lmol=lbb(l)
               amult=1.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=0.5d0
               endif
               if(ifno) then
                  totalbb=totalbb+amult*(pbb(l)*(sij(i,ik)*sij(j,jl)+
     .                   sij(j,ik)*sij(i,jl))-dsqrt(pbb(l))*(sij(i,il)*
     .                   sij(j,jk)+sij(j,il)*sij(i,jk)))
               else
                  totalbb=totalbb+amult*pbb(l)*(sij(i,ik)*sij(j,jl)+
     .                   sij(j,ik)*sij(i,jl)-sij(i,il)*sij(j,jk)-
     .                   sij(j,il)*sij(i,jk))
               endif
            enddo
            bbij(i,j)=totalbb
            covb(i,j)=totalbb-popb(i)*popb(j)
         enddo
         covb(i,i)=covb(i,i)+popb(i)
         stdev(i)=dsqrt(covb(i,i))
      enddo
      write(*,'(//,'' beta spin covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                                      (covb(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                           (covb(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      else
      do i=1,nbas
         do j=1,nbas
            covb(i,j)=cova(i,j)
            bbij(i,j)=aaij(i,j)
         enddo
      enddo
      endif
      do i=1,nbas
         do j=1,nbas
            totalab=0.0d0
            do l=1,nab
               imol=iab(l)
               jmol=jab(l)
               kmol=kab(l)
               lmol=lab(l)
               amult=2.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               if(uhf.and.ifci) then
                  jl=indij(jmol+n_alpha)+lmol+n_alpha
               endif
               if(ik.eq.jl) then
                  amult=1.0d0
               endif
               if(uhf) then
                  amult=0.5d0*amult
               endif
               totalab=totalab+amult*pab(l)*(sij(i,ik)*sij(j,jl)+
     .                                       sij(j,ik)*sij(i,jl))
            enddo
            covab(i,j)=aaij(i,j)+bbij(i,j)+totalab-rho_bas(i)*rho_bas(j)
            abij(i,j)=totalab
         enddo
         covab(i,i)=covab(i,i)+rho_bas(i)
         stdev(i)=dsqrt(covab(i,i))
      enddo
      write(*,'(//)')
      write(*,'(//,'' total covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                                      (covab(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nbas) then
            j2=nbas
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,nbas
            write(*,'(1x,a16,10f6.2)')name_basin(i),
     .                          (covab(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      write(*,'(/,2x,''elapsed time'',f9.2,''sec.'')') elapsed
      endif
      if(rho_pop) then
      total_charge=total
      do i=1,natoms
         totala=0.0d0
         totalb=0.0d0
         totalaa=0.0d0
         totalbb=0.0d0
         totalab=0.0d0
         do l=1,na
            lmol=ia(l)
            totala=totala+pa(l)*tij(i,lmol)
         enddo
         do l=1,naa
            imol=iaa(l)
            jmol=jaa(l)
            kmol=kaa(l)
            lmol=laa(l)
            amult=2.0d0
            ik=indij(imol)+kmol
            if(imol.gt.kmol) then
               ik=indij(kmol)+imol
            endif
            jl=indij(jmol)+lmol
            if(jmol.gt.lmol) then
               jl=indij(lmol)+jmol
            endif
            il=indij(imol)+lmol
            if(imol.gt.lmol) then
               il=indij(lmol)+imol
            endif
            jk=indij(jmol)+kmol
            if(jmol.gt.kmol) then
               jk=indij(kmol)+jmol
            endif
            if(ik.eq.jl) then
               amult=1.0d0
            endif
            if(ifno) then
               totalaa=totalaa+amult*(paa(l)*tij(i,ik)*tij(i,jl)-
     .                        dsqrt(paa(l))* tij(i,il)*tij(i,jk))
            else
               totalaa=totalaa+amult*paa(l)*(tij(i,ik)*tij(i,jl)-
     .                                    tij(i,il)*tij(i,jk))
            endif
         enddo
         if(nb.ne.0) then
            do l=1,nb
               lmol=ib(l)
               totalb=totalb+pb(l)*tij(i,lmol)
            enddo
            totalbb=0.0d0
            do l=1,nbb
               imol=ibb(l)
               jmol=jbb(l)
               kmol=kbb(l)
               lmol=lbb(l)
               amult=2.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=1.0d0
               endif
               if(ifno) then
                  totalbb=totalbb+amult*(pbb(l)*tij(i,ik)*tij(i,jl)-
     .                            dsqrt(pbb(l))*tij(i,il)*tij(i,jk))
               else
                  totalbb=totalbb+amult*pbb(l)*(tij(i,ik)*tij(i,jl)-
     .                                    tij(i,il)*tij(i,jk))
               endif
            enddo
         else
            totalb=totala
            totalbb=totalaa
         endif
*     sum for opposite spin
         do l=1,nab
            imol=iab(l)
            jmol=jab(l)
            kmol=kab(l)
            lmol=lab(l)
            amult=2.0d0
            ik=indij(imol)+kmol
            if(imol.gt.kmol) then
               ik=indij(kmol)+imol
            endif
            jl=indij(jmol)+lmol
            if(jmol.gt.lmol) then
               jl=indij(lmol)+jmol
            endif
            if(uhf.and.ifci) then
               jl=indij(jmol+n_alpha)+lmol+n_alpha
            endif
            if(ik.eq.jl) then
               amult=1.0d0
            endif
            totalab=totalab+amult*pab(l)*tij(i,ik)*tij(i,jl)
         enddo
         total=totala+totalb
         rho_at(i)=total
         popa(i)=totala
         popb(i)=totalb
         if(uhf) then
            popab(i)=totalab
         else
            popab(i)=2.0d0*totalab
         endif
         popaa(i)=totalaa
         popbb(i)=totalbb
         arg=popaa(i)+popbb(i)+popab(i)+rho_at(i)*(1.0d0-rho_at(i))
         bij(i,i)=arg
         stdev(i)=dsqrt(arg)
         if(rho_at(i).eq.0.0d0) then
           rho_at(i)=1.0d-3
         endif
      enddo
      total=0.0d0
      write(*,'(//,t3,''AIM  ANALYSIS '')')
      write(*,'(//,t3,''basin'',t23,''vol.'',t31,''pop.'',t39,"pab",
     .          t47,"paa",t54,"pbb",t62,''sigma2'',
     .          t69,''std. dev.'',//)')
      do i=1,natoms
         write(*,'(i3,1x,a16,f6.2,6(2x,f6.2))') i,name_atom(i),
     .                 vol_at(i),rho_at(i),popab(i),popaa(i),popbb(i),
     .                 bij(i,i),stdev(i)
         total=total+rho_bas(i)
      enddo
      total=0.0d0
      do i=1,natoms
         total=total+rho_at(i)
      enddo
      write(*,'(//,'' sum of atomic populations '',f12.6,/)') total
      if(.not.ifci) then
      write(*,'(//,'' orbital contributions '',/)')
      j1=-9
      j2=0
      jlimit=.true.
      imol=nmol
      if(fukui) imol=nmol-1
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.nmol) then
            j2=nmol
            jlimit=.false.
         endif
         write(*,'(//,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
c        allocate(pop_orb(natoms))
         do i=1,natoms
            if(j2.le.imol) then
               do j=j1,j2
                  ij=nmol*(j-1)+(3*j-j*j)/2
                  pop_orb(j)=(occ(1,j)+occ(2,j))*tij(i,ij)
               enddo
            else
               ij=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
               pop_orb(nmol)=tij(i,ij)
            endif
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                                (pop_orb(j),j=j1,j2)
        enddo
      enddo
      if(fukui) then
         if(open_shell) then
         write(*,'(//,'' condensed Fukui indexes '',/)')
         j1=nmol*(imol-1)+(3*imol-imol*imol)/2
         j2=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
         j3=nmol*(n_beta-1)+(3*n_beta-n_beta*n_beta)/2
         j4=nmol*n_beta+(n_beta-n_beta**2+2)/2
         if(multiplicity.gt.2) then
            write(*,'(t22,''f-'',t28,''f+'',t34,"fs+",t39,"fs-",//)')
            do i=1,nbas
            write(*,'(1x,a16,4f6.2)')name_atom(i),tij(i,j1),tij(i,j4),
     .                                0.5d0*(tij(i,j2)+tij(i,j3)),
     .                                0.5d0*(tij(i,j1)+tij(i,j4))
            enddo
            else
            write(*,'(t22,''f-'',t28,''f+'',t34,"fs+",t39,//)')
            do i=1,nbas
            write(*,'(1x,a16,3f6.2)')name_atom(i),tij(i,j1),tij(i,j4),
     .                                0.5d0*(tij(i,j2)+tij(i,j3))
            enddo
         endif
         else
         write(*,'(//,'' condensed Fukui indexes '',/)')
         j1=nmol*(imol-1)+(3*imol-imol*imol)/2
         j2=nmol*(nmol-1)+(3*nmol-nmol*nmol)/2
         write(*,'(t22,''f-'',t28,''f+'',t34,''f0'',t39,''Df'',//)')
         do i=1,natoms
            write(*,'(1x,a16,4f6.2)')name_atom(i),tij(i,j1),tij(i,j2),
     .                                0.5d0*(tij(i,j1)+tij(i,j2)),
     .                                  tij(i,j1)-tij(i,j2)
         enddo
        endif
      endif
      endif
      if(open_shell) then
         write(*,'(//,t1,''basin'',t20,''integrated spin density'',//)'
     .            )
         totals=0.0d0
         do i=1,natoms
            total=0.5d0*(popa(i)-popb(i))
            write(*,'(i3,1x,a20,2x,f12.6)') i,name_atom(i),total
            totals=totals+total
         enddo
         write(*,'(//,'' sum of spin densities '',f12.6)') totals
      endif
      do i=1,natoms
         sum=0.0d0
         do j=1,natoms
            totalaa=0.0d0
            do l=1,naa
               imol=iaa(l)
               jmol=jaa(l)
               kmol=kaa(l)
               lmol=laa(l)
               amult=1.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=0.5d0
               endif
               if(ifno) then
                  totalaa=totalaa+amult*(paa(l)*(tij(i,ik)*tij(j,jl)+
     .                   tij(j,ik)*tij(i,jl))-dsqrt(paa(l))*(tij(i,il)*
     .                   tij(j,jk)+tij(j,il)*tij(i,jk)))
               else
                  totalaa=totalaa+amult*paa(l)*(tij(i,ik)*tij(j,jl)+
     .                   tij(j,ik)*tij(i,jl)-tij(i,il)*tij(j,jk)-
     .                   tij(j,il)*tij(i,jk))
               endif
            enddo
            aaij(i,j)=totalaa
            cova(i,j)=aaij(i,j)-popa(i)*popa(j)
         enddo
         cova(i,i)=aaij(i,i)+popa(i)*(1.0d0-popa(i))
         stdev(i)=dsqrt(cova(i,i))
      enddo
      write(*,'(//,'' alpha spin covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                                      (cova(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                           (cova(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      if(open_shell) then
      do i=1,natoms
         do j=1,natoms
            totalbb=0.0d0
            do l=1,nbb
               imol=ibb(l)
               jmol=jbb(l)
               kmol=kbb(l)
               lmol=lbb(l)
               amult=1.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               il=indij(imol)+lmol
               if(imol.gt.lmol) then
                  il=indij(lmol)+imol
               endif
               jk=indij(jmol)+kmol
               if(jmol.gt.kmol) then
                  jk=indij(kmol)+jmol
               endif
               if(ik.eq.jl) then
                  amult=0.5d0
               endif
               if(ifno) then
                  totalbb=totalbb+amult*(pbb(l)*(tij(i,ik)*tij(j,jl)+
     .                   tij(j,ik)*tij(i,jl))+dsqrt(pbb(l))*(tij(i,il)*
     .                   tij(j,jk)+tij(j,il)*tij(i,jk)))
               else
                  totalbb=totalbb+amult*pbb(l)*(tij(i,ik)*tij(j,jl)+
     .                   tij(j,ik)*tij(i,jl)-tij(i,il)*tij(j,jk)-
     .                   tij(j,il)*tij(i,jk))
               endif
            enddo
            bbij(i,j)=totalbb
            covb(i,j)=bbij(i,j)-popb(i)*popb(j)
         enddo
         covb(i,i)=covb(i,i)+popb(i)
         stdev(i)=dsqrt(covb(i,i))
      enddo
      write(*,'(//,'' beta spin covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                                      (covb(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                           (covb(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      else
      do i=1,natoms
         do j=1,natoms
            bbij(i,j)=aaij(i,j)
         enddo
      enddo
      endif
      do i=1,natoms
         do j=1,natoms
            totalab=0.0d0
            do l=1,nab
               imol=iab(l)
               jmol=jab(l)
               kmol=kab(l)
               lmol=lab(l)
               amult=2.0d0
               ik=indij(imol)+kmol
               if(imol.gt.kmol) then
                  ik=indij(kmol)+imol
               endif
               jl=indij(jmol)+lmol
               if(jmol.gt.lmol) then
                  jl=indij(lmol)+jmol
               endif
               if(uhf.and.ifci) then
                   jl=indij(jmol+n_alpha)+lmol+n_alpha
               endif
               if(ik.eq.jl) then
                   amult=1.0d0
               endif
               if(uhf) then
                  amult=0.5d0*amult
               endif
                totalab=totalab+amult*pab(l)*(tij(i,ik)*tij(j,jl)+
     .                                        tij(j,ik)*tij(i,jl))
            enddo
            abij(i,j)=totalab
            covab(i,j)=aaij(i,j)+bbij(i,j)+abij(i,j)-rho_at(i)*rho_at(j)
         enddo
         covab(i,i)=covab(i,i)+rho_at(i)
         stdev(i)=dsqrt(covab(i,i))
      enddo
      write(*,'(//)')
      write(*,'(//,'' total covariance matrix'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                                      (covab(i,j),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      write(*,'(//,'' correlation coefficients'')')
      j1=-9
      j2=0
      jlimit=.true.
      do while (jlimit)
         j1=j1+10
         j2=j2+10
         if(j2.ge.natoms) then
            j2=natoms
            jlimit=.false.
         endif
         write(*,'(/,16x,10i6)')(j,j=j1,j2)
         write(*,'(/)')
         do i=1,natoms
            write(*,'(1x,a16,10f6.2)')name_atom(i),
     .                          (covab(i,j)/(stdev(i)*stdev(j)),j=j1,j2)
         enddo
      enddo
      write(*,'(//)')
      call xtime(t)
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      write(*,'(/,2x,''elapsed time'',f9.2,''sec.'')') elapsed
      endif
      end
