       function coulomb(x1j,x2j,x1j1,x2j1,ldim)
       implicit none
c regularized coulomb action with gaussian integration
       include 'pimc.par'
       integer il,ldim,l,iu
       real*8 x1j(ldim),x2j(ldim),x1j1(ldim),x2j1(ldim),coulomb,d,dl
       real*8 sg,erf,arg,erfc,bessi0,dd4,dd
       external erf,erfc,bessi0,bessi1
       include 'pimc.cm'

       coulomb=0.d0
       do iu=1,iord
        d=0.d0
        do l=1,ldim
         dl=pbc(
     &          x1j(l)*up2(iu)+x1j1(l)*um2(iu),
     &          x2j(l)*up2(iu)+x2j1(l)*um2(iu),l
     &         )
         d=d+dl*dl
        enddo
        d=sqrt(d)
c       if (d.lt.cutr) then
         sg=sigma2*dsqrt(1.d0-4.d0*xnode(iu)**2)
         arg=d/sg
         coulomb=coulomb+wg(iu)*erf(arg)/d
!       endif
       enddo
       coulomb=coulomb*charge*tau
       return
       end

       function dcoulomb(x1j,x2j,x1j1,x2j1,ldim)
       implicit none
c regularized coulomb action with gaussian integration
       include 'pimc.par'
       integer il,ldim,l,iu
       real*8 x1j(ldim),x2j(ldim),x1j1(ldim),x2j1(ldim),dcoulomb,d,dl
       real*8 sg,erf,erfc,bessi0,bessi1,dd4,dd,dd2
       external erf,erfc,bessi0,bessi1
       include 'pimc.cm'

       dcoulomb=0.d0
       do iu=1,iord
        d=0.d0
        do l=1,ldim
         dl=pbc(
     &          x1j(l)*up2(iu)+x1j1(l)*um2(iu),
     &          x2j(l)*up2(iu)+x2j1(l)*um2(iu),l
     &         )
         d=d+dl*dl
        enddo
        dd2=d
        sg=sigma2*sigma2*(1.d0-4.d0*xnode(iu)**2)
        d=d/sg
        sg=dsqrt(sg)*spi
        dcoulomb=dcoulomb-wg(iu)*dexp(-d)/sg
       enddo
       dcoulomb=dcoulomb*charge
       return
       end
 


       function coulvir(x1j,x2j,x1j1,x2j1,x11,x12,ldim)
       implicit none
c regularized coulomb action with gaussian integration
       include 'pimc.par'
       integer il,ldim,l,iu
       real*8 x1j(ldim),x2j(ldim),x1j1(ldim),x2j1(ldim),coulvir,d
       real*8 sg,dl(mdim),dl2,ajab,bjab,erfc,x11(ldim),x12(ldim)
       real*8 arg,dla,dd2,bessi0,bessi1,dd,dd4
       external erfc,bessi0,bessi1
       include 'pimc.cm'

       coulvir=0.d0
       do iu=1,iord
        dd2=0.d0
        ajab=0.d0
        do l=1,ldim
         dl(l)=pbc(
     &         x1j(l)*up2(iu)+x1j1(l)*um2(iu),
     &         x2j(l)*up2(iu)+x2j1(l)*um2(iu),l
     &            )
         dd2=dd2+dl(l)*dl(l)
         ajab=ajab+pbc(x11(l),x12(l),l)*dl(l)
        enddo
        if(ndim.eq.3)then
         if (dd2.lt.cutr**2) then
          sg=sigma2*sigma2*(1.d0-4.d0*xnode(iu)**2)
          dl2=dd2/sg
          dla=alp**2*dd2
          sg=dsqrt(sg)
          d=dsqrt(dd2)
          arg=d/sg
          bjab=2.d0/spi*(dexp(-dl2)/sg-alp*dexp(-dla))
     &        +(erfc(arg)-erfc(alp*d))/d
          ajab=ajab/dd2
          coulvir=coulvir+wg(iu)*bjab*(1.d0-ajab)
         endif
        else
         if (dd2.lt.cutr**2) then
	  sg=sigma2*sigma2*(1.d0-4.d0*xnode(iu)**2)
	  dd4=dd2/sg/2.d0
	  dla=alp**2*dd2
	  sg=dsqrt(2.d0*sg)
	  bjab= 2.d0*sq2pi*(dd4/sg*(bessi1(dd4)-bessi0(dd4))
     &         -alp**3*dd2*(bessi1(dla)-bessi0(dla)))
	  ajab=ajab/dd2
	  coulvir=coulvir+wg(iu)*bjab*(1.d0-ajab)
	 endif
        endif
       enddo
       return
       end
