subroutine cvirtab
c construct the virial estimator vectors A_j^\alpha B_j^{\alpha}
c for the coulomb interaction term
c 6/11/01 

      implicit none
      include 'pimc.par'
      integer js,js1,l,iu,n,m
      real*8 dl(mdim),dd2,d,sg2,sg,dl2,arg,dla,bjab,ddd
      real*8 bessi0,bessi1,erfc
      external erfc,bessi0,bessi1
      include 'pimc.cm'

      do n=1,np
       do js=1,nslices
        js1=iwrap(js-1)
        do l=1,ndim
         aja(l,n,js)=0.d0
         bja(l,n,js)=0.d0
        enddo
        do m=1,n-1
         do iu=1,iord
c computes L_i^{nm}(u)=dl in PBC
          dd2=0.0
          do l=1,ndim
           dl(l)=pbc(x(l,n,js1)*um2(iu)+x(l,n,js)*up2(iu),
     &               x(l,m,js1)*um2(iu)+x(l,m,js)*up2(iu),l)
           dd2=dd2+dl(l)**2
          enddo
          if (dd2.lt.cutr2) then
           dla=alp**2*dd2
           sg2=sigma2*sigma2*(1.d0-4.d0*xnode(iu)**2)
           dl2=dd2/sg2
           sg=dsqrt(sg2)
           d=dsqrt(dd2)
           arg=d/sg
           if (ndim.eq.3) then
            bjab=(2.d0/spi*(dexp(-dl2)/sg-alp*dexp(-dla))
     &          +(erfc(arg)-erfc(alp*d))/d)/dd2
           elseif(ndim.eq.2) then
            dl2=dl2/2.d0
            sg=sq2*sg
            bjab= (2.d0*sq2pi*(dl2/sg*(bessi1(dl2)-bessi0(dl2)
     &          -alp**3*dd2*(bessi1(dla)-bessi0(dla)))))/dd2
           endif
           do l=1,ndim
            ddd=wg(iu)*bjab*dl(l)
            aja(l,n,js)=aja(l,n,js)+ddd
            aja(l,m,js)=aja(l,m,js)-ddd
            bja(l,n,js)=bja(l,n,js)+ddd*xnode(iu)
            bja(l,m,js)=bja(l,m,js)-ddd*xnode(iu)
           enddo
          endif
         enddo   ! over iu
        enddo    ! over m
       enddo     ! over js
      enddo      ! over n
      return
      end