********************************************************************************
*                                                                              *
*  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 ### focksl ###
c
c     Calculates the exchange potentials as the local Slater approximation.
c
      subroutine focksl(iorb,psi,pot,excp,e,f0,f1,f2,f4,wk0,wk1,wk2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'	
      dimension psi(*),pot(*),excp(*),e(*),f0(*),f1(*),f2(*),f4(*),
     &	    wk0(*),wk1(*),wk2(*)	

      iborb=i1b(iorb)
      ngorb=i1si(iorb)
      ibpot=i2b(iorb)
      ngpot=i2si(iorb)	
      norb2=norb*norb			

      if (idbg(97).ne.0) then
         write(*,*) 'iorb,iborb,ngorb,ibpot,ngpot'
         write(*,*)  iorb,iborb,ngorb,ibpot,ngpot
      endif

      ipe=mgx(6,iorb)
      ipex=ipe-2*(ipe/2)
      
      if (ipe.eq.0) then
	 call dcopy (mxsize,f0,1,wk0,1)
	 call daxpy (mxsize,eng(iorb),f1,1,wk0,1)
      else

c        e enter the expression with minus sign which is already
c        incorporated in e

         w=dble(ipe*ipe)
         call dcopy (mxsize,f0,1,wk0,1)
         call daxpy (mxsize,eng(iorb),f1,1,wk0,1)
         call daxpy (mxsize,w,e,1,wk0,1)
      endif

      if(idbg(97).ne.0) then
         write(*,*) 'eng(iorb)',eng(iorb)
         write(*,*) 'focksl: f0'
         call pmtx (nni,nmu(1),f0,1,1,incrni,incrmu)
         write(*,*) 'focksl: f1'
         call pmtx (nni,nmu(1),f1,1,1,incrni,incrmu)
         write(*,*) 'focksl: f2'
         call pmtx (nni,nmu(1),f2,1,1,incrni,incrmu)
         write(*,*) 'focksl: f4'
         call pmtx (nni,nmu(1),f4,1,1,incrni,incrmu)
         write(*,*) 'fock: wk0'
         call pmtx (nni,nmu(1),wk0,1,1,incrni,incrmu)
      endif

      do i=1,mxsize
         wk1 (i)=0.d0
         wk2 (i)=0.d0
         foc2(i)=0.d0
         excp(i)=0.d0
      enddo	

      if (nel.gt.1) then
         if (norb.eq.1) then

c     contribution from coulomb potential of one sigma orbital

c     24/04/01  
c     Valentin Karasiev spotted and fixed an error in handling 
c     two-electron cases (for a single orbital systems)
c           call dcopy (ngpot,pot(ibpot),1,wk1,1)
c     12/12/02 jkob added
c           call prod (mxsize,f2,wk1)
c     in order to treat systems like LiH+2 properly
      
            oc=occ(iorb)
            call daxpy (ngpot,oc,pot(ibpot),1,wk1,1)
            call prod (mxsize,f2,wk1)
         else
            
c           add contributions from coulomb and off-diagonal lagrange mult.

            do iorb1=1,norb
               iborb1=i1b(iorb1)	   	 
               ngorb1=i1si(iorb1)	 
               ibpot1=i2b(iorb1)
               ngpot1=i2si(iorb1)

               ideng=iorb1+norb*(iorb-1)
               coo=occ(iorb1)
               
c              in the local exchange approx. the coulomb potential includes
c              also the contribution from the orbital

               if (idbg(95).ne.0) then
                  write(*,*) 'focksl: iorb1,ibpot1,coo',iorb1,ibpot1,coo
               endif
               call daxpy (ngpot1,coo,pot(ibpot1),1,wk1,1)

               if (iorb.ne.iorb1.and.engo(ideng).ne.0.d0) then
                  do i=1,ngorb1
                     foc2(i)=foc2(i)+engo(ideng)*f4(i)*psi(iborb1+i-1)
                  enddo		
               endif
            enddo

c        contributions from the local slater exchange

            do iorb1=1,norb
               iborb1=i1b(iorb1)	   	 
               ngorb1=i1si(iorb1)	 	 
               coo=occ(iorb1)
               do i=1,ngorb1
                  wk2(i)=wk2(i)+coo*psi(iborb1+i-1)*psi(iborb1+i-1)
               enddo
            enddo

c           multiply the exchange potentia by f4 to make it commensurate 
c           with the coulomb potntial

            do i=1,mxsize
               wk2(i)=-3.d0/2.d0*alphaf
     &                *(3.d0/pii*wk2(i))**(1.d0/3.d0)*f4(i)
            enddo		
            
            if(idbg(97).ne.0) then
               write(*,*) 'vfock: alphaf,sl1,sl2 ',alphaf,sl1,sl2
               write(*,*) 'nni,nmu,incrni,incrmu',
     &              nni,nmu(1),incrni,incrmu
               write(*,*) 'vfock: slater potential'
               call pmtx (nni,nmu(1),wk2,1,1,incrni,incrmu)
               
               write(*,*) 'focksl: coulomb potential'
               call pmtx (nni,nmu(1),wk1,1,1,incrni,incrmu)
               
               write(*,*) 'focksl: one-electron contr.'
               call pmtx (nni,nmu(1),wk0,1,1,incrni,incrmu)
            endif
            
c           store the local exchange pot. in excp array
            
            call dcopy (ngorb1,wk2,1,excp(1),1)
            
c           add the local exchange to the coulomb potential
            
            call add (mxsize,wk2,wk1)
            
c           multiply coulomb/exchange potentials and off-diagonal Lagrange
c           multipliers by f2
         
            call prod (mxsize,f2,wk1)
            call prod (mxsize,f2,foc2)
         endif
c        nel.gt.1 
      endif

c     add one electron and coulomb pot. contributions

      call add (mxsize,wk0,wk1)
      call dcopy (mxsize,wk1,1,foc1,1)

      return
      end