subroutine coulinit(inow)
      implicit none
! initialize potential and action for coulomb potential
      include 'tas.cm'
      real*8 tpiell(3),wtk(0:mnshex),b(20),work(2000),cusp(6)
     .,u(mnc,mnc),vmad,eps,pot,dpot,x,xcut
     .,totcharge,eh,const,arg,th,q,s,tn,vk,ukk
      integer l,iexp(6),ifcon(6),it,jt,mpoly,nderv,nshex,i,k
     .,inow,ifl
      character p(1)*8

! mh
      real*8 ulradd(mnc,mnc)
      real*8 kfm(2),a3,s0,ss,y
      real*8 yq,uqadd,wy,wu
      real*8 dx,dx2,dxn,ddpot,rr
      integer nk1,isp
! end mh
      save ifl
      include 'cpbc.cm'
      data ifl/0/
      ifl=ifl+1


      nshex=mnshex
      call fillk(rknorm(0,inow),wtk,kmult(0,inow),nkact(inow),nshex
     .    ,ndim,volume(inow))

!     if(root)write (6,*)' nelects =',nelects,' volume',volume(inow)
!    .,' inow=',inow,' nkact(inow)=',nkact(inow)
      do k=1,nshex
        call murpa(ntypes,rknorm(k,inow),charge,hbs2m,u,mnc,rho,root)
        do it=1,ntypes
        do jt=1,ntypes
          ulr(k,it,jt,2,inow)=u(it,jt)/nelects
c ulr(k,1) is the coulomb potential
          ulr(k,it,jt,1,inow)=(charge(it)*charge(jt)*pi*2*(ndim-1))
     .     /(volume(inow)*rknorm(k,inow)**(ndim-1))
       enddo
       enddo
!      if(debug.and.root)
!    .  write (6,*)rknorm(k,inow),
!    .   ((ulr(k,it,jt,2,inow),it=1,ntypes),jt=1,ntypes)
      enddo

! mh
      
      kfm(1)=0.d0
      kfm(2)=0.d0
      if(nspins(ite).gt.2) then
           if(root)write(6,*)'coulinit nspins(e)=',nspins(ite)
           if(root)write(6,*)'is bigger than 2 !'
           call parallel_end()
      end if
      do isp=1,nspins(ite)
        a3=3.d0*volume(inow)/(4.d0*pi*nppss(isp,ite))
        kfm(isp)=(9.d0*pi/(2.d0*a3))**(1.d0/ndim)
      end do
      do k=1,nshex
       s0=0.d0  ! the ideal gas structure factor
       do isp=1,nspins(ite)
        if(nppss(isp,ite).gt.0) then
         y=0.5d0*(rknorm(k,inow))/kfm(isp)
         if(y.lt.1.d0) then
          if(ndim.eq.3) then
           ss=0.5*y*(3.d0-y*y)
          elseif(ndim.eq.2) then
           ss=(2.d0/pi)*(dasin(y)+y*dsqrt(1.d0-y**2))
          endif
         else
          ss=1.d0
         endif
         s0=s0+nppss(isp,ite)*ss/nelects
        endif
       enddo
       do it=1,ntypes
         do jt=1,ntypes
           call uanalytique(k,kfm,s0,it,jt,uqadd,yq,wu,wy,inow)
           ulr(k,it,jt,3,inow)=-yq ! backflow
           ulradd(it,jt)=uqadd ! self/extra jastrow
           ulr(k,it,jt,4,inow)=wy ! 3body (backflow)
           ulr(k,it,jt,5,inow)=wu ! part of 3-body potential (jastrow)
         enddo
       enddo
       do it=1,ntypes
         do jt=1,ntypes
           ulr(k,it,jt,2,inow)=ulr(k,it,jt,2,inow)
     .                        +ulradd(it,jt)
!           if(root.and.debug)write(50+it*ntypes+jt+inow*10,'(6f15.6)')
!     .            (ulr(k,it,jt,l,inow),l=2,5)
         end do
       end do
      end do

! end mh

      totcharge=0.d0
      do it=1,ntypes
      totcharge=totcharge+charge(it)*ncomps(it)
      enddo
      totcharge=totcharge/charge(1)
!     if(root)write (6,*)' total charge ',totcharge

c Find optimal division into k and r space.
      do l=1,5 ! loop over the potential, jastrow and bf factor
       do it=1,ntypes
        do jt=1,ntypes
         xcut=cutr(inow)
         mpoly=10            !  order of the polynomial
         nderv=2            ! number of derivatives zeroed at cutr
         vmad=0.d0
         nk1=nkact(inow)+1
         eps=charge(it)*charge(jt)
         ulr(0,it,jt,l,inow)=0.d0 !  set k=0 term to zero
         if(l.eq.1) then    ! potential
          iexp(1)=0
          cusp(1)=eps
          ifcon(1)=1
          if(ndim.eq.3)vmad=-2.837297479d0*eps/ell(1)
          if(ndim.eq.2)vmad=-3.90026492d0*eps/ell(1)
          if(pname(it).eq.'p'.and.pname(jt).eq.'p') then
           if (npslices.gt.1.and.ifew.gt.0) xcut=cutrp
          endif
         elseif(l.eq.2) then  ! jastrow
          iexp(2)=1
          eh=hbs2m(it)+hbs2m(jt)
          if(eh.gt.0.d0) then
             cusp(2)=-eps/((ndim-1)*eh)
          else
             cusp(2)=0.d0
          endif
          ifcon(2)=1
          vmad=0.d0
! mh
         elseif(l.eq.3) then ! backflow
          nderv=3
          iexp(3)=1 ! constant term at origin
          cusp(3)=0.d0 ! cusp and ifcon set the first derivative at origin to zero
          ifcon(3)=1
          vmad=0.d0
!          nk1=nbact+1
          nk1=min(nbact,nkact(inow))+1
        elseif(l.eq.4) then ! backflow-part of 3-body
          iexp(4)=1 ! constant term at origin
          cusp(4)=0.d0 ! cusp and ifcon set the first derivative at origin to zero
          ifcon(4)=1
          vmad=0.d0
!          nk1=n3act+1  
          nk1=min(n3act,nkact(inow))+1
          nderv=3 ! since we want the derivative
        elseif(l.eq.5) then   ! jastrow-part of 3-body
          iexp(5)=1
          cusp(5)=0.d0
          ifcon(5)=1
          vmad=0.d0
!          nk1=n3act+1
          nk1=min(n3act,nkact(inow))+1
          nderv=3
        end if
!        if(root.and.debug) write(6,*)'========================'
!        if(root.and.debug) write(6,*)'l=',l,' it=',it,' jt=',jt
        call fitpn(mpoly,ulr(0,it,jt,l,inow),rknorm(0,inow),wtk
     .   ,nshex+1,nk1,cusp(l),b,work,iexp(l),volume(inow)
     +   ,xcut,ndim,ifcon(l),nderv,vmad,root)

! mh end
         if(l.eq.1) then
!           if(root.and.debug)
!     .         write (6,*) 'const=',it,jt,ulr(0,it,jt,l,inow)
           if(it+jt.eq.3) const=ulr(0,it,jt,l,inow)
           if(it+jt.eq.2)ulr(0,1,1,1,inow)=
     .      0.5d0*ulr(0,1,1,1,inow)*totcharge**2
           if(it.eq.jt)ulr(0,1,1,1,inow)=ulr(0,1,1,1,inow)
     .       +vmad*ncomps(it)*.5d0
         endif
c now evaluate potential
         do k=0,lptable
           rr=dtable(inow)*k
           if(k.eq.0) rr=dtable(inow)*.001d0
           x=rr/xcut
           pot=0.d0
           dpot=0.d0
           ddpot=0.d0
           if(x.le.1.d0) then
              do i=mpoly-nderv,1,-1
                 pot=b(i)+x*pot  ! function
              enddo
              do i=mpoly-nderv,2,-1
                 dpot=(i-1)*b(i)+x*dpot  ! derivative
              enddo
              do i=mpoly-nderv,3,-1
                 ddpot=(i-1)*(i-2)*b(i)+x*ddpot  ! derivative
              enddo
              dx=x-1.d0
              dx2=dx*dx
              dxn=1.d0
              if(nderv.ne.2) dxn=dx**(nderv-2)
              ddpot=(dx2*ddpot+2*nderv*dx*dpot+nderv*(nderv-1)*pot)*dxn
              dpot=(nderv*pot+dx*dpot)*dx*dxn
              pot=pot*dx2*dxn ! mulitply by dx**nderv
           endif
c     for Coulomb potential fit usr*r
           if(l.eq.1) then
              usr(k,1,it,jt,l,inow)=pot*xcut
              usr(k,2,it,jt,l,inow)=dpot*dtable(inow)
! write out the sr potential table
!             if(it.eq.1.and.jt.eq.1) then
!              if(ifl.eq.1) then
!               write (93,'(3g18.8)') rr,pot/x
!              endif
!             endif
           elseif(l.eq.2)then
              urpa(k,it,jt,1)=pot
              urpa(k,it,jt,2)=dpot/xcut ! dtable will be added in ptrans
           elseif(l.eq.3) then ! Backflow
              ubackflow(k,it,jt,1)=dpot/(rr*xcut)
              ubackflow(k,it,jt,2)=
     .           (ddpot/xcut-dpot/rr)/(rr*xcut)
           elseif(l.eq.4) then ! Backflow for 3-body
              u3back(k,it,jt,1)=dpot/(rr*xcut)
              u3back(k,it,jt,2)=
     .           (ddpot/xcut-dpot/rr)/(rr*xcut)
           elseif(l.eq.5) then ! Jastrow for 3-body
              u3jas(k,it,jt,1)=dpot/(rr*xcut)
              u3jas(k,it,jt,2)=
     .           (ddpot/xcut-dpot/rr)/(rr*xcut)
           end if
         enddo
        end do
       enddo
      enddo
      do it=1,ntypes
        do jt=1,ntypes
          call spline(usr(0,1,it,jt,1,inow),lptable+1,lptable+1,0) ! potential
        end do
      end do
! correlation function will be done in ptrans
      pot_expon=1 ! means potential is table divided by r.
      vtail=vimage+ulr(0,1,1,1,inow)
!     if(root)write(6,*)'coulinit: potential tail = ',vtail
      call ptrans(0,0,p,inow)

!      if(root.and.debug)
!     .   write(6,*)'coulinit: potential tail = ',ulr(0,1,1,1,inow)
!      call ptrans(0,0,p,inow)
!     do l=1,2
!     do k=1,nkact
!       if(root)write (6,'(4f14.6)')((ulr(k,it,jt,l),it=1,ntypes),jt=1,ntypes)
!     enddo
!     enddo
!     ulr(0,1,1,1)=const  !DEBUG for ken
!
! END-POINT correction to the pair action for PIMC
!
      if(npslices.gt.1.and.ifew.gt.0) then
!      taup=31.579
!      hbar2m=2.7231e-4
       do k=1,nshex
        do it=1,ntypes
         if(pname(it).eq.'p') then
          vk=ulr(k,it,it,1,inow) ! long range part of Coulomb pot.
          ulr(k,it,it,1,inow)=0.d0
          ulr(k,it,it,2,inow)=0.d0
         endif
        enddo
        arg=1.d0+2.d0*vk*dble(nprotons)/(hbar2m*rknorm(k,inow)**2)
        th=0.5d0*rknorm(k,inow)**2*hbar2m*taup
        if(arg.le.0.d0)then
         q=dsqrt(-arg)
         if(q*th.ge.0.5d0*pi) then
          if(root) write (6,*) 'coulinit: q*th.ge.pi/2'
          call parallel_end()
         end if
         tn=dtan(th*q)
         s=q*(1.d0-q*tn)/(q+tn)
        elseif(arg.eq.0.d0)then
         s=1.d0/(1.d0+th)
        elseif(arg.gt.0.d0)then
         q=dsqrt(arg)
         tn=dtanh(th*q)
         s=q*(1.d0+q*tn)/(q+tn)
        endif
        ukk=(-1.d0+s)/dble(nprotons)
        ukep(k,1,inow)=vk
        ukep(k,2,inow)=ukk
        ukep(k,3,inow)=-rknorm(k,inow)**2*hbar2m*ukk
     .                 *(1.d0+0.5d0*dble(nprotons)*ukk)
!       write(94,'(4g15.8)')rknorm(k,inow),vk,ukk,ukep(k,3,inow)
       enddo
      endif

      end