********************************************************************************
*                                                                              * 
*  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 ### tfpot ###  
c
c     
      subroutine tfpot (psi,pot,excp,f2,f4,wk0)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),f2(*),f4(*),wk0(*)

      if (imethod.eq.2.or.ini.eq.4) return
      
      if (ini.ne.6) then
         
c        Initialization of Coulomb potentials
         
         print *,'... initializing Coulomb potentials ...'
         
c        loop over orbitals
         
         slim=1.0d0
         zc1= 2.0d0
         zc2= 2.0d0
         do iorb=1,norb
            ngrid=i1si(iorb)	 	 
            ishift=i1b(iorb)-1	
            ez1=eza1(iorb)
            ez2=eza2(iorb)	  	
            
            ch1=-1.0d0
            ch2=-1.0d0
            if (ez1.lt.precis) ch1=-2.0d0
            if (ez2.lt.precis) ch2=-2.0d0
            
            crt1=abs(co1(iorb))
            crt2=abs(co2(iorb))
            
c           loop over grid points

            do imu=1,mxnmu
               inioff=(imu-1)*nni
               vxit=vxi(imu)
               do in=1,nni
                  igp=ishift+inioff+in
                  vetat=veta(in)
                  if (imu.eq.1.and.in.eq.1) then
                     pot(igp)=0.d0
                  else
                     pot(igp)=thofe(r,vetat,vxit,zc1,ch1,slim)*crt1
     &                    +thofe(r,(-vetat),vxit,zc2,ch2,slim)*crt2
                  endif
               enddo
            enddo
            call prod(ngrid,f4,pot(i1b(iorb)))               
         enddo
      endif

c     Initialization of exchange potentials
      
      do i=1,length3
         excp(i)=0.d0
      enddo

      if (imethod.eq.3) then
         call slaterp(psi,pot,excp,f2,f4,wk0)
         print *,'... initializing Slater exchange potential ...'
         return
      elseif (iform.eq.1.or.iform.eq.3) then
         
         print *,'... initializing exchange potentials ...'

c        Initialization of exchange potentials that are all kept in memory
         
c        loop over orbitals
         
         do iorb1=1,norb
            do iorb2=iorb1,norb
               k=iorb1+iorb2*(iorb2-1)/2
               if (iorb1.eq.iorb2.and.ll(iorb2).eq.0) goto 50
               ngrid=i3si(k)
               ishift=i3b(k)-1
               
c              loop over grid points
           
               do imu=1,mxnmu
                  if (imu.eq.1) goto 40
                  inioff=(imu-1)*nni
                  vxit=vxi(imu)
                  do in=1,nni
                     if (in.eq.1.or.in.eq.nni) goto 30
                     igp=ishift+inioff+in
                     vetat=veta(in)
                     ra1=(r/2.d0)*(vxit+vetat)
                     ra2=(r/2.d0)*(vxit-vetat)

                     if (ra1.gt.1.d-8) then
                        excp(igp)=co1(iorb1)/ra1
                     else
                        excp(igp)=0.d0
                     endif
                     if (ra2.gt.1.d-8) then
                        excp(igp)=excp(igp)+co2(iorb1)/ra2
                     endif
 30                  continue
                  enddo
 40               continue	
               enddo
               call prod(ngrid,f4,excp(i3b(k)))               
               if (iorb1.eq.iorb2) goto 50
               if (ll(iorb1).eq.0.or.ll(iorb2).eq.0) goto 50
               ishift=ishift+ngrid
               call dcopy(ngrid,excp(i3b(k)),1,excp(i3b(k)+ngrid),1)
 50            continue
            enddo
         enddo

      elseif (iform.eq.0.or.iform.eq.2) then

         print *,'... initializing exchange potentials ...'

         call prepwexch
         do iorb1=1,norb
            do iorb2=1,i3nexcp(iorb1)
               k=i3breck(iorb1,iorb2)
               i3beg=i3brec(iorb1,iorb2)
               iorb2t=i3orb2(k)
               i3b(k)=i3beg        
               ngrid=i3si(k)
               if (ilc(k).eq.1) then
                  irec=i3xpair(iorb1,iorb2)
                  ishift=i3b(k)-1
                  
c                 loop over grid points
                  
                  if (iwexch(k).ne.0) then
                     do imu=1,mxnmu
                        if (imu.eq.1) goto 41
                        inioff=(imu-1)*nni
                        vxit=vxi(imu)
                        do in=1,nni
                           if (in.eq.1.or.in.eq.nni) goto 31
                           igp=ishift+inioff+in
                           vetat=veta(in)
                           ra1=(r/2.d0)*(vxit+vetat)
                           ra2=(r/2.d0)*(vxit-vetat)
                          
                           if (ra1.gt.1.d-8) then
                              excp(igp)=co1(iorb1)/ra1
                           else
                              excp(igp)=0.d0
                           endif
                           if (ra2.gt.1.d-8) then
                              excp(igp)=excp(igp)+co2(iorb1)/ra2
                           endif
 31                        continue
                        enddo
 41                     continue	
                     enddo
                  endif
                  call prod(ngrid,f4,excp(i3b(k)))               
               elseif (ilc(k).eq.2) then
                  irec=i3xpair(iorb1,iorb2)
                  call dcopy(ngrid,excp(i3beg-ngrid),1,
     &                 excp(i3beg),1)
               endif
            enddo
            
c           save exchange potential involving iorb1 on disk 
            call wtdexch(iorb1,excp)
         enddo
      endif
      return
      end