********************************************************************************
*                                                                              *
*  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 ### multi ###
c
c     Calculates the value for the multipole expansion for the exchange 
c     potential in practical infinity.

      subroutine multi(i,j,mt,ipc,pe)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension dome(10)

      rr=sqrt(vxisq(j)+vetasq(i)-1.d0)		
      costh=0.d0
      if (abs(rr).gt.precis) then
         costh=veta(i)*vxi(j)/rr	 
      endif

      costh2=costh*costh
      costh3=costh2*costh
      costh4=costh2*costh2
      rr=rr*r2
      pe=0.d0

      if (mt.eq.0) then
c        m=0
         dome( 1)=costh
         dome( 2)=( 3.d0*costh*dome(1) - 1.d0)	  / 2.d0
         dome( 3)=( 5.d0*costh*dome(2) - 2.d0*dome(1))/ 3.d0
         dome( 4)=( 7.d0*costh*dome(3) - 3.d0*dome(2))/ 4.d0
         dome( 5)=( 9.d0*costh*dome(4) - 4.d0*dome(3))/ 5.d0
         dome( 6)=(11.d0*costh*dome(5) - 5.d0*dome(4))/ 6.d0
         dome( 7)=(13.d0*costh*dome(6) - 6.d0*dome(5))/ 7.d0
         dome( 8)=(15.d0*costh*dome(7) - 7.d0*dome(6))/ 8.d0

c        dome( 9)=(17.d0*costh*dome(8) - 8.d0*dome(7))/ 9.d0
c        dome(10)=(19.d0*costh*dome(9) - 9.d0*dome(8))/10.d0

         pe=      dome(8)*exc8 (ipc) /rr
         pe=(pe + dome(7)*exc7 (ipc))/rr
         pe=(pe + dome(6)*exc6 (ipc))/rr

         pe=(pe + dome(5)*exc5 (ipc))/rr
         pe=(pe + dome(4)*exche(ipc))/rr
         pe=(pe + dome(3)*excoc(ipc))/rr
         pe=(pe + dome(2)*excqu(ipc))/rr
         pe=(pe + dome(1)*excdi(ipc))/rr/rr
         return	

c     pe=(3.5d01*costh4-3.d01*costh2+3.d0)*exche(ipc)*calp(1)/rr
c     pe=(pe+(5.d0*costh3-3.d0*costh)*excoc(ipc)*calp(2))/rr
c     pe=(pe+(3.d0*costh2-1.d0)*excqu(ipc)*calp(2))/rr
c     rx2=rr*rr
c     pe=(pe+costh*excdi(ipc))/rx2

      elseif (mt.eq.1) then
c        m=1
         sini2=abs(1.d0-costh2)
         sini=sqrt(sini2)
         pe=sini*(7.d0*costh3-3.d0*costh)*exche(ipc)*calp(3)/rr
         pe=(pe+sini*(5.d0*costh2-1.d0)*excoc(ipc)*calp(4))/rr
         pe=(pe+sini*costh*excqu(ipc)*calp(5))/rr
         rx2=rr*rr
         pe=(pe+sini*excdi(ipc)*calp(12))/rx2
         
         return
         
      elseif (mt.eq.2) then
c        m=2
         sini2=1.d0-costh2
         pe=sini2*(7.d0*costh2-1.d0)*exche(ipc)*calp(6)/rr
         pe=(pe+sini2*costh*excoc(ipc)*calp(7))/rr
         rx3=rr*rr*rr
         pe=(pe+sini2*excqu(ipc)*calp(8))/rx3
         return
      elseif (mt.eq.3) then
c        m=3
         sini2=abs(1.d0-costh2)
         sini=sqrt(sini2)
         sini3=sini2*sini
         pe=sini3*costh*exche(ipc)*calp(9)/rr
         rx4=rr*rr*rr*rr
         pe=(pe+sini3*excoc(ipc)*calp(10))/rx4
         return
      elseif (mt.eq.4) then
c        m=4
         sini2=1.d0-costh2
         sini4=sini2*sini2
         rx5=rr*rr*rr*rr*rr
         pe=sini4*exche(ipc)*calp(11)/rx5
      endif
      return
      end