********************************************************************************
*                                                                              * 
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### summary ####
c
      subroutine summary (iscf,cw_orb,cw_coul,cw_exch,
     &     cw_suppl,cw_sctch)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension cw_orb(*),cw_coul(*),cw_exch(*),
     &     cw_suppl(*),cw_sctch(*)

c     write functions to disk file if nobckup < 10000

      if (nobckup.ge.0) then
         call wtdisk(cw_orb,cw_coul,cw_exch)
      endif

c     calculate final total energy
      
      write (*,*)
      call toten (cw_orb,cw_coul,cw_exch,                         
     &     cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b(13)),
     &     cw_suppl(i4b(14)),
     &     cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),
     &     cw_sctch(i5b( 4)))

      write(iout6,6110)  etot
      write(iout6,6100)  evt
      if (islat.eq.0) write(iout6,6120)  virrat
      etot_final=etot


c     write final orbital energies and normalization factors
      
      write(*,*)
      write(iout6,1000)
      do i=1,norb
         write(iout6,1002) iorn(i),bond(i),gut(i),eng(i),area(i)
      enddo  
      write(*,*)

c     check how an error in the orbital normalization factors influence
c     the accuracy of the total energy

      do iorb=1,norb
         ibeg = i1b (iorb)
         ngrid= i1si(iorb)
         xnorm=area(iorb)
         zarea = sqrt(xnorm)
         call dscal (ngrid,zarea,cw_orb(ibeg),1)
      enddo  

c     calculate final total energy
      
      write (*,*)
      call toten (cw_orb,cw_coul,cw_exch,                         
     &     cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b(13)),
     &     cw_suppl(i4b(14)),
     &     cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),
     &     cw_sctch(i5b( 4)))

      write(iout6,6111)  abs(etot_final-etot),
     &     abs(etot_final/etot-1.d0)*100

      do iorb=1,norb
         ibeg = i1b (iorb)
         ngrid= i1si(iorb)
         xnorm=area(iorb)
         zarea = 1.d0/sqrt(xnorm)
         call dscal (ngrid,zarea,cw_orb(ibeg),1)
      enddo  
 
      
c     check orthohonality

      if (idbg(210).ne.0) then
         write(*,*) 'checking orthogonality:'
         iprt=2
         do iorb=1,norb
            i=norb+1-iorb
            call ortho (i,cw_orb,
     &           cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &           cw_sctch(i5b(1)),iprt)
         enddo
      endif

c     check symmetry of orbitals in the homonuclear case
      
      if (abs(z1-z2).lt.homolevl) call homo_check(cw_orb)

c      call orbtails (cw_orb,cw_coul,cw_exch,cw_suppl,cw_sctch)
      
      if (idbg(495).ne.0) then
         write(*,*)
         write(*,*) 'Euclidean norms of (T+V(n)+V-E)|i>:'
         do iorb=1,norb
            call locenergy (iorb,cw_orb,cw_coul,cw_exch,
     &           cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b( 9)),
     &           cw_suppl(i4b(13)),cw_suppl(i4b(14)),
     &           cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),
     &           cw_sctch(i5b( 4)))
         enddo
      endif

c     check multipole expansion
      if (idbg(302).ne.0) then
         call chasympot (cw_coul,cw_exch)
      endif

c     calculate multipole moments and expectation values

      call propet (iscf,cw_orb,cw_coul,cw_exch,cw_suppl(i4b( 9)),
     &     cw_suppl(i4b(14)),
     &     cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)),          
     &     cw_sctch(i5b( 4)),cw_sctch(i5b( 5)),cw_sctch(i5b( 6)))
c     
      write(*,'(//,"////////////////////////////////////////////",
     &     "////////////////////////////////////////////////////",
     &     "////////////////")')

c     cpu time statistics
      
      tcpu=trelax+tortho+trayl+tmomen

      write(*,'(/," cpu summary (sec):")')
      write(*,5001) trelax
      write(*,5002) tortho
      write(*,5003) trayl
      write(*,5004) tmomen
      write(*,5005) tcpu
      if (iscf.ne.0) write(*,5006) tcpu/dble(iscf)
      
 1000 format(19x,'  orbital energies       normalization factors')  
 1002 format(1x,i3,1x,a8,1x,a1,2d25.13)
 5001 format(4x,'relaxation of orb. and pot.        ', f12.2)
 5002 format(4x,'normalization+orthogonalization    ', f12.2)
 5003 format(4x,'lagrange multipliers               ', f12.2)
 5004 format(4x,'multipole moments                  ', f12.2)
 5005 format(4x,'sum total                          ', f12.2)
 5006 format(4x,'duration of single scf iteration   ', f12.2)
 6100 format(1x,'total electronic energy: ',d25.13)
 6110 format(1x,'total energy:            ',d25.13)
 6111 format(1x,'total energy uncertaininty due to orbital norms',
     &     ' not being equal 1:'/,
     &     '  absolute   +/-',d8.2,/,
     &     '  relative   +/-',d8.2,'%')
 6120 format(1x,'virial ratio:            ',d25.13)

      return
      end