********************************************************************************
*                                                                              * 
*  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 ### totensl ###
c
c     Calculates total energy for the local exchange approximation
c
      subroutine totensl (psi,pot,excp,e,f0,wgt1,wgt2,wk0,wk1,wk2,wk3)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'	
      dimension psi(*),pot(*),excp(*),e(*),f0(*),wgt1(*),wgt2(*),
     &     wk0(*),wk1(*),wk2(*),wk3(*)  

      data eps /1.d-7/

      engt(1)=eng(1)

      if (nel.lt.2) return

      engt(1)=0.d0
      
c     calculate first contributions from one particle operators and
c     Coulomb potential contributions within the same shell

      if (idbg(70).ne.0) then
	 write(*,*) 'orb        t        v(n)      t+v(n)     ',
     &        ' v(c-within shell)'
      endif

      norb2=norb*norb
      woneel=0.d0
      wdcoul=0.d0

      do iorb=1,norb
         iborb=i1b(iorb)
         isiorb=i1si(iorb)
         nmut=i1mu(iorb)
         ibpot=i2b(iorb)	  
         isipot=i2si(iorb)
         oc=occ(iorb)
         
         ipe=mgx(6,iorb)
         ipex=ipe-2*(ipe/2)

c        calculate derivatives over mu and ni variables by means of matrix
c        multiplication

         if (ipex.eq.0) then
	    isym= 1
         else
	    isym=-1
         endif
         
         call fill (nni,nmut,isym,psi(iborb),wk3)
         call difni (nmut,wk3,wk0,wk1,wk2)
         call refill(nni,nmut,wk1,wk0)
         call difmu (nmut,wk3,wk2)
         call refill(nni,nmut,wk0,wk2)

c        add derivatives over mu and ni

         call add (isiorb,wk0,wk1)

c        add contribution from phi part of laplasian
c        e enters the expression with minus sign which is already
c        incorporated in e

         if (ipe.ne.0) then
	    w=dble(ipe*ipe)
	    call dcopy (isiorb,e,1,wk0,1)
	    if (ipe.ne.1) then
               call dscal (isiorb,w,wk0,1)
	    endif
	    call prod  (isiorb,psi(iborb),wk0)
	    call add   (isiorb,wk0,wk1)
         endif
         
         call dcopy (isiorb,wk1,1,wk2,1)
         call prod  (isiorb,psi(iborb),wk2)
         w=ddot(isiorb,wgt1,1,wk2,1)
         vk(iorb)=w
         
         call dcopy (isiorb,f0,1,wk0,1)
         call prod  (isiorb,psi(iborb),wk0)
         
         call prod2 (isiorb,psi(iborb),wk0,wk2)
         w=ddot(isiorb,wgt1,1,wk2,1)
         vn(iorb)=w

         call add	 (isiorb,wk0,wk1)
         call prod (isiorb,psi(iborb),wk1)
         
         w = ddot(isiorb,wgt1,1,wk1,1)
         woneel=woneel+oc*w
      enddo

c     contribution from coulomb interaction within the same shell

      if (isipot.ne.isiorb) then
         write(iout6,*) 'coulomb potentials and orbitals have to ',
     &        'defined on the same number of subgrids'
         stop 'toten'
      endif				

c     calculate the coulomb potential contribution from all the orbitals
c     (include 1/2 factor )

      do i=1,mxsize
         wk2(i)=0.d0
      enddo
      
      do iorb=1,norb
         ibpot=i2b(iorb)
         isipot=i2si(iorb)	
         oc=occ(iorb)/2.d0
         call daxpy (isipot,oc,pot(ibpot),1,wk2,1)
      enddo

c     contribution from coulomb interaction

      wndc =0.d0
      do iorb=1,norb
         iborb=i1b(iorb)
         isiorb=i1si(iorb)	  
         call prod2 (isiorb,psi(iborb),psi(iborb),wk0)
         call prod (isiorb,wk2,wk0)
         call dscal (isiorb,occ(iorb),wk0,1)
         w=ddot(isiorb,wgt2,1,wk0,1)
         wndc=wndc+w
      enddo
      
c     contribution from local exchange interaction

      do i=1,mxsize
         wk1(i)=0.d0
         wk2(i)=0.d0
      enddo	  
      
      do iorb=1,norb
         iborb=i1b(iorb)
         isiorb=i1si(iorb)	  
         call exocc (iorb,ocup,ocdown)
         call prodas (isiorb,ocup ,psi(iborb),psi(iborb),wk1)
         call prodas (isiorb,ocdown,psi(iborb),psi(iborb),wk2)
      enddo
      
      do i=1,mxsize
         wk1(i)=wk1(i)**(4.d0/3.d0)
         wk2(i)=wk2(i)**(4.d0/3.d0)
      enddo	  

      call add (mxsize,wk1,wk2)
      
c     since the local exchange potentai is formed from the orbitals directly
c     multiply it by f4 to get rid of 2/(r*cosh(xi)) already present in
c     wgt2

      xrr=r/2.d0
      do i=1,mxnmu
         ii=(i-1)*nni
         do j=1,nni
	    wk0(ii+j)=xrr*vxi(i)
         enddo
      enddo

      call prod (mxsize,wk0,wk2)

      w=ddot(mxsize,wgt2,1,wk2,1)
      wex=-9.d0/4.d0*alphaf*(3/(4.d0*pii))**(1.d0/3.d0) * w
      
      evt=woneel+wndc+wex
      engt(1)=evt
      etot=evt+z1*z2/r

 6100 format(1x,'total electronic energy: ',d25.13)
 6110 format(1x,'total            energy: ',d25.13)

      if (idbg(70).ne.0) then
         write(*,*)'totesl - woneel,wndc,wex,etotal ',
     &        woneel,wndc,wex,evt
      endif

      return
      end