********************************************************************************
*                                                                              *
*  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 ### fock ###
c
c     Calculates the Fock potential for a given orbital.
c
      subroutine fock (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		
      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(*,*) 'fock: f0'
         call pmtx (nni,nmu(1),f0,1,1,incrni,incrmu)
         write(*,*) 'fock: f1'
         call pmtx (nni,nmu(1),f1,1,1,incrni,incrmu)
         write(*,*) 'fock: f2'
         call pmtx (nni,nmu(1),f2,1,1,incrni,incrmu)
         write(*,*) 'fock: 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
      enddo	

      if (nel.gt.1) then
         if (norb.eq.1) then
            
c           contribution from coulomb potential of one sigma orbital

            call dcopy (ngpot,pot(ibpot),1,wk1,1)
         else

c           add contributions from coulomb and exchange potentials

            do iorb1=1,norb
               iborb1=i1b(iorb1)
               ngorb1=i1si(iorb1)
               kex=iorb+norb*(iorb1-1)
               ideng=iorb1+norb*(iorb-1)
               
               if (iorb.le.iorb1) then
                  idexp=iorb+iorb1*(iorb1-1)/2
               else
                  idexp=iorb1+iorb*(iorb-1)/2
               endif
               
               ibpot1=i2b(iorb1)
               ngpot1=i2si(iorb1)
               ibexp=i3b(idexp)
               ngexp=i3si(idexp)
               
               if (idbg(95).ne.0) then
                  write(*,*) 'focksl: iorb1,ibpot1,coo',iorb1,ibpot1,coo
               endif
               
               coo=occ(iorb1)
               if (iorb.eq.iorb1) coo=coo-1.d0
               
               call daxpy (ngpot1,coo,pot(ibpot1),1,wk1,1)
               
               if (iorb1.ne.iorb)  then
                  coo=gec(kex)
                  call dcopy (ngexp,excp(ibexp),1,wk2,1)
                  call dscal (ngexp,coo,wk2,1)
                  
                  if (ilc(idexp).gt.1) then
                     coo=gec(kex+norb2)
                     call daxpy (ngexp,coo,excp(ibexp+ngexp),1,wk2,1)
                  endif
                  
                  if (engo(ideng).ne.0.d0) then
                     call daxpy (ngexp,engo(ideng),f4,1,wk2,1)
                  endif
                  
                  call prod (ngorb1,psi(ibpot1),wk2)
                  call add  (ngorb1,wk2,foc2)
               else
                  if ((ipe.gt.0).and.(ilc(idexp).gt.0)) then
                     coo=gec(kex)
                     call dcopy (ngexp,excp(ibexp),1,wk2,1)
                     ngrid=max(ngexp,ngorb1)
                     call prod  (ngorb1,psi(iborb1),wk2)
                     call dscal (ngrid,coo,wk2,1)
                     call add (ngrid,wk2,foc2)
                  endif
               endif
            enddo
         endif
         
         if (idbg(97).ne.0) then
            write(*,*) 'fock: coulomb potential'
            call pmtx (nni,nmu(1),wk1,1,1,incrni,incrmu)
            write(*,*) 'fock: one-electron contr.'
            call pmtx (nni,nmu(1),wk0,1,1,incrni,incrmu)
c           write(*,*) ' fock0 - f4   '
c           call pmtx (nni,nmi,f4   ,incrni,incrmi)
         endif
         
c        multiply coulomb and exchange pot. by f2
         
         call prod (mxsize,f2,wk1)
         call prod (mxsize,f2,foc2)
      endif

c  add one electron and coulomb pot. contributions
c  are coulomb and exchange pot zero at infinity?

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

      return
      end