subroutine lrange(rnew,gradu,u,pe,del2u,task,jslice,inow) ! adds k-space part to v u and pe
      implicit none
       include 'tas.cm'
      real*8  rnew(ndim,nparts),gradu(ndim,nparts),srhok
     .,u,pe,del2u,cons,rhok(mnkv2,mnc)
      integer j,ks,l,i,nk,k,task,it,jt,j0,inow,ispec,jslice,ip,ie,ifind
!     real*8 rhokpp(mnkv2),pw(mnkv2,mpart)
! fixes pwm and rhok

! mh
      integer i1,m,i2,nk1
      real*8 xdummy,xi,ehb,xa2,xd1,xd2
     .,xa,xi1,xi2
     

      real*8 rkmatqee(ndim,ndim),rkmatqep(ndim,ndim)
     .,rkvqee(ndim),rkvqep(ndim),rk2ee,rk2ep,xaee,xaep,xiee,xiep
      real*8 rkvw1ee(ndim),rkvw1ep(ndim),rkvw2ee(ndim),rkvw2ep(ndim)
     .,rkmatw1ee(ndim,ndim),rkmatw2ee(ndim,ndim)
     .,rkmatw1ep(ndim,ndim),rkmatw2ep(ndim,ndim)

      real*8 h1,h2,force1,force2,dforce1,dforce2
     .,v,w,d2v,d2w

      common/c3body/h1(ndim,ndim,mpart,mpart)
     .,h2(ndim,ndim,mpart,mpart)
     .,v(ndim,mpart,mpart),w(ndim,mpart,mpart)
     .,d2v(ndim,mpart,mpart),d2w(ndim,mpart,mpart)
     .,force1(ndim,mpart),dforce1(ndim,mpart)
     .,force2(ndim,mpart),dforce2(ndim,mpart)


      ip=ifind('p',pname,ntypes)
c     if(npslices.eq.1) then
! also protons for classical 
       call CompRhok(ip,rnew(1,nfpty(ip)),rhokp(1,jslice,inow)
     .              ,pwmatp(1,1,jslice,inow),inow)
c     endif
      ie=ifind('e',pname,ntypes)
      call CompRhok(ie,rnew(1,nfpty(ie)),rhoke(1,inow),pwmate,inow) 

! set rhok and pwm for local calculation
      do k=1,nvects2(inow)
       rhok(k,ie)=rhoke(k,inow)
       rhok(k,ip)=rhokp(k,jslice,inow)
!      if(rhokp(k,jslice,inow).ne.rhokpp(k)) 
!    .  write(*,*)inow,nvects2(inow),k,rhokp(k,jslice,inow),rhokpp(k)
       do i=1,ncomps(ie)
        j=i-1+nfpty(ie)
        pwmat(k,j)=pwmate(k,i)
       enddo
       do i=1,ncomps(ip)
        j=i-1+nfpty(ip)
        pwmat(k,j)=pwmatp(k,i,jslice,inow)
!       if(pwmatp(k,i,jslice,inow).ne.pw(k,i))
!    .   write(*,*)inow,ncomps(ip),k,i,pwmatp(k,i,jslice,inow),pw(k,i)
       enddo
      enddo
       
      j0=1
      do ks=1,nkact(inow)  ! long-range potential
       nk=kmult(ks,inow)-kmult(ks-1,inow)
       do it=1,ntypes
        do jt=1,it
         ispec=0
         if(pname(it).eq.'p'.and.pname(jt).eq.'p') ispec=1
! uncomment the following line for unbiased PIMC of proton
!        if(npslices.gt.1.and.ispec.eq.1) goto 10 ! skip for p-p with PIMC
         srhok=0.d0
         j=j0
         do k=kmult(ks-1,inow)+1,kmult(ks,inow)
          srhok=srhok+rhok(j,it)*rhok(j,jt)+rhok(j+1,it)*rhok(j+1,jt)
          j=j+2
         enddo
         if(it.ne.jt) srhok=2*srhok
         if(ifpair.ne.0)u=u+srhok*ulr(ks,it,jt,2,inow)
         pe=pe+srhok*ulr(ks,it,jt,1,inow) 
         if(task.gt.0.and.ifpair.ne.0) then
          if(it.eq.jt)srhok=srhok-nk*ncomps(it)
          del2u=del2u-(hbs2m(it)+hbs2m(jt))*rknorm(ks,inow)**2
     .                                     *ulr(ks,it,jt,2,inow)*srhok
          j=j0
          do k=kmult(ks-1,inow)+1,kmult(ks,inow)
           do l=1,ndim ! for the gradient
            cons=2.d0*ulr(ks,it,jt,2,inow)*rkcomp(l,k,inow)
            do i=nfpty(it),nlpty(it)
            gradu(l,i)=gradu(l,i)-cons*(pwmat(j+1,i)*rhok(j,jt)
     .                                -pwmat(j,i)*rhok(j+1,jt))
            enddo
            if(it.ne.jt) then
             do i=nfpty(jt),nlpty(jt)
             gradu(l,i)=gradu(l,i)-cons*(pwmat(j+1,i)*rhok(j,it)
     .                                -pwmat(j,i)*rhok(j+1,it))
             enddo
            endif
           end do
! mh end

              j=j+2
            enddo
         endif
10       continue
        enddo
       enddo
       j0=j0+2*nk
      enddo
!   mh
!      nk1=nbact
      nk1=min(nbact,nkact(inow))
      j0=1
      do ks=1,nk1
        rk2ep=rknorm(ks,inow)**2*hbs2m(ite)
        rk2ee=rk2ep*2.d0
         
        j=j0
        do k=kmult(ks-1,inow)+1,kmult(ks,inow)
          do l=1,ndim
           xdummy=2.d0*rkcomp(l,k,inow) 
           rkvqee(l)=xdummy*ulr(ks,ite,ite,3,inow)
           rkvqep(l)=xdummy*ulr(ks,ite,itp,3,inow)
           do m=1,l
             xdummy=2.d0*rkcomp(l,k,inow)*rkcomp(m,k,inow)
             rkmatqee(l,m)=xdummy*ulr(ks,ite,ite,3,inow)
             rkmatqep(l,m)=xdummy*ulr(ks,ite,itp,3,inow)
           end do
          end do
          do i1=nfpty(ite),nlpty(ite)
             xiee= pwmat(j,i1)*rhok(j+1,ite)
     .          -pwmat(j+1,i1)*rhok(j,ite)
             xaee=(1.d0-pwmat(j+1,i1)*rhok(j+1,ite)
     .               -pwmat(j,i1)*rhok(j,ite))
     
             xiep= pwmat(j,i1)*rhok(j+1,itp)
     .            -pwmat(j+1,i1)*rhok(j,itp)
             xaep=-(pwmat(j+1,i1)*rhok(j+1,itp)
     .             +pwmat(j,i1)*rhok(j,itp))

             do l=1,ndim
                xi1=rkvqee(l)*xiee
                xi2=rkvqep(l)*xiep
                qp(l,i1)=qp(l,i1)+xi1+xi2
                bmat(l,i1)=bmat(l,i1)-xi1*rk2ee
     .                               -xi2*rk2ep
                do m=1,l
                   amat(l,m,i1,i1)=amat(l,m,i1,i1)
     .                  +rkmatqee(l,m)*xaee            
     .                  +rkmatqep(l,m)*xaep
                end do
             end do

!             if(ks.le.nbact) then
               do i2=nfpty(ite),i1-1
                xa= pwmat(j,i2)*pwmat(j,i1)
     .             +pwmat(j+1,i2)*pwmat(j+1,i1)

                do l=1,ndim
                   do m=1,l
                      amat(l,m,i1,i2)=amat(l,m,i1,i2)
     .                   +rkmatqee(l,m)*xa
                   end do
                end do
               end do
!             end if
          end do
          j=j+2
        end do
        j0=j0+2*(kmult(ks,inow)-kmult(ks-1,inow))
      end do
      if(ifbacka.gt.0) then
        do i1=nfpty(ite),nlpty(ite)
         do i2=nfpty(ite),i1-1
            do l=1,ndim
               amat(l,l,i2,i1)=amat(l,l,i1,i2)
               do m=1,l-1
                  amat(m,l,i1,i2)=amat(l,m,i1,i2)
                  amat(l,m,i2,i1)=amat(l,m,i1,i2)
                  amat(m,l,i2,i1)=amat(l,m,i1,i2)
               end do
            end do
         end do
        end do
      end if              


      if(ifbackflow.ne.0) then
      do i=1,nparts      ! symmetrize amat along diagonal
        do l=2,ndim
        do m=1,l-1
           amat(m,l,i,i)=amat(l,m,i,i)
        enddo
        enddo
      enddo
      endif


!   mh
      j0=1
      nk1=min(n3act,nkact(inow))
      do ks=1,nk1
        rk2ep=rknorm(ks,inow)**2*hbs2m(ie)
        rk2ee=rk2ep*2.d0
         
        j=j0
        do k=kmult(ks-1,inow)+1,kmult(ks,inow)
          do l=1,ndim
           xdummy=2.d0*rkcomp(l,k,inow) 
           rkvw1ee(l)=xdummy*ulr(ks,ite,ite,4,inow)
           rkvw2ee(l)=xdummy*ulr(ks,ite,ite,5,inow)
           rkvw1ep(l)=xdummy*ulr(ks,ite,itp,4,inow)
           rkvw2ep(l)=xdummy*ulr(ks,ite,itp,5,inow)
           do m=1,l
             xdummy=2.d0*rkcomp(l,k,inow)*rkcomp(m,k,inow)
             rkmatw1ee(l,m)=xdummy*ulr(ks,ite,ite,4,inow)
             rkmatw1ep(l,m)=xdummy*ulr(ks,ite,itp,4,inow)
             rkmatw2ee(l,m)=xdummy*ulr(ks,ite,ite,5,inow)
             rkmatw2ep(l,m)=xdummy*ulr(ks,ite,itp,5,inow)
           end do
          end do
          do i1=nfpty(ie),nlpty(ie)
             xiee= pwmat(j,i1)*rhok(j+1,ie)
     .          -pwmat(j+1,i1)*rhok(j,ie)
             xaee=(1.d0-pwmat(j+1,i1)*rhok(j+1,ie)
     .               -pwmat(j,i1)*rhok(j,ie))

             xiep= pwmat(j,i1)*rhok(j+1,ip)
     .            -pwmat(j+1,i1)*rhok(j,ip)
             xaep=-(pwmat(j+1,i1)*rhok(j+1,ip)
     .             +pwmat(j,i1)*rhok(j,ip))

                do l=1,ndim
                   xi1=rkvw1ee(l)*xiee
                   xi2=rkvw1ep(l)*xiep
                   force1(l,i1)=force1(l,i1)+xi1+xi2
                   dforce1(l,i1)=dforce1(l,i1)-xi1*rk2ee
     .                                        -xi2*rk2ep
                   xi1=rkvw2ee(l)*xiee
                   xi2=rkvw2ep(l)*xiep
                   force2(l,i1)=force2(l,i1)+xi1+xi2
                   dforce2(l,i1)=dforce2(l,i1)-xi1*rk2ee
     .                                        -xi2*rk2ep
                   do m=1,l
                      h1(l,m,i1,i1)=h1(l,m,i1,i1)
     .                     +rkmatw1ee(l,m)*xaee
     .                     +rkmatw1ep(l,m)*xaep
                      h2(l,m,i1,i1)=h2(l,m,i1,i1)
     .                     +rkmatw2ee(l,m)*xaee
     .                     +rkmatw2ep(l,m)*xaep
                   end do
                end do
                do i2=nfpty(ie),i1-1
                   xi= pwmat(j,i1)*pwmat(j+1,i2)
     .                -pwmat(j+1,i1)*pwmat(j,i2)
                   do l=1,ndim
                      xi1=rkvw1ee(l)*xi
                      v(l,i1,i2)=v(l,i1,i2)+xi1
                      v(l,i2,i1)=v(l,i2,i1)-xi1
                      xi1=xi1*rk2ee
                      d2v(l,i1,i2)=d2v(l,i1,i2)-xi1
                      d2v(l,i2,i1)=d2v(l,i2,i1)+xi1
                      xi1=rkvw2ee(l)*xi
                      w(l,i1,i2)=w(l,i1,i2)+xi1
                      w(l,i2,i1)=w(l,i2,i1)-xi1
                      xi1=xi1*rk2ee
                      d2w(l,i1,i2)=d2w(l,i1,i2)-xi1
                      d2w(l,i2,i1)=d2w(l,i2,i1)+xi1
                   end do
                   xa= pwmat(j,i2)*pwmat(j,i1)
     .                +pwmat(j+1,i2)*pwmat(j+1,i1)
                   do l=1,ndim
                      do m=1,l
                         h1(l,m,i1,i2)=h1(l,m,i1,i2)
     .                     +rkmatw1ee(l,m)*xa
                         h2(l,m,i1,i2)=h2(l,m,i1,i2)
     .                     +rkmatw2ee(l,m)*xa
                      end do
                   end do
                end do
                do i2=nfpty(ip),nlpty(ip)
                   xi= pwmat(j,i1)*pwmat(j+1,i2)
     .                -pwmat(j+1,i1)*pwmat(j,i2)
                   do l=1,ndim
                      xi1=rkvw1ep(l)*xi
                      v(l,i1,i2)=v(l,i1,i2)+xi1
                      v(l,i2,i1)=v(l,i2,i1)-xi1
                      xi1=xi1*rk2ep
                      d2v(l,i1,i2)=d2v(l,i1,i2)-xi1
                      d2v(l,i2,i1)=d2v(l,i2,i1)+xi1
                      xi1=rkvw2ep(l)*xi
                      w(l,i1,i2)=w(l,i1,i2)+xi1
                      w(l,i2,i1)=w(l,i2,i1)-xi1
                      xi1=xi1*rk2ep
                      d2w(l,i1,i2)=d2w(l,i1,i2)-xi1
                      d2w(l,i2,i1)=d2w(l,i2,i1)+xi1
                   end do
                   xa= pwmat(j,i2)*pwmat(j,i1)
     .                +pwmat(j+1,i2)*pwmat(j+1,i1)
                   do l=1,ndim
                      do m=1,l
                        h1(l,m,i1,i2)=h1(l,m,i1,i2)
     .                   +rkmatw1ep(l,m)*xa
                        h2(l,m,i1,i2)=h2(l,m,i1,i2)
     .                   +rkmatw2ep(l,m)*xa
                      end do
                   end do
                end do
          end do
             do i1=nfpty(ip),nlpty(ip)
                xiep= pwmat(j,i1)*rhok(j+1,ie)
     .            -pwmat(j+1,i1)*rhok(j,ie)
                xaep=-(pwmat(j+1,i1)*rhok(j+1,ie)
     .             +pwmat(j,i1)*rhok(j,ie))
                do l=1,ndim
                   xi=rkvw1ep(l)*xiep
                   force1(l,i1)=force1(l,i1)+xi
                   dforce1(l,i1)=dforce1(l,i1)-xi*rk2ep
                   xi=rkvw2ep(l)*xiep
                   force2(l,i1)=force2(l,i1)+xi
                   dforce2(l,i1)=dforce2(l,i1)-xi*rk2ep
                   do m=1,l
                      h1(l,m,i1,i1)=h1(l,m,i1,i1)
     .                +rkmatw1ep(l,m)*xaep
                      h2(l,m,i1,i1)=h2(l,m,i1,i1)
     .                +rkmatw2ep(l,m)*xaep
                   end do
                end do
             end do

          j=j+2
        end do
        j0=j0+2*nk1
      end do

      if(if3bodya.gt.0) then
        do i1=nfpty(ie),nlpty(ie)
         do i2=nfpty(ie),i1-1
            do l=1,ndim
               h1(l,l,i2,i1)=h1(l,l,i1,i2)
               h2(l,l,i2,i1)=h2(l,l,i1,i2)
               do m=1,l-1
                  h1(m,l,i1,i2)=h1(l,m,i1,i2)
                  h1(l,m,i2,i1)=h1(l,m,i1,i2)
                  h1(m,l,i2,i1)=h1(l,m,i1,i2)
                  h2(m,l,i1,i2)=h2(l,m,i1,i2)
                  h2(l,m,i2,i1)=h2(l,m,i1,i2)
                  h2(m,l,i2,i1)=h2(l,m,i1,i2)
               end do
            end do
         end do
         do i2=nfpty(ip),nlpty(ip)
            do l=1,ndim
               h1(l,l,i2,i1)=h1(l,l,i1,i2)
               h2(l,l,i2,i1)=h2(l,l,i1,i2)
               do m=1,l-1
                  h1(m,l,i1,i2)=h1(l,m,i1,i2)
                  h1(l,m,i2,i1)=h1(l,m,i1,i2)
                  h1(m,l,i2,i1)=h1(l,m,i1,i2)
                  h2(m,l,i1,i2)=h2(l,m,i1,i2)
                  h2(l,m,i2,i1)=h2(l,m,i1,i2)
                  h2(m,l,i2,i1)=h2(l,m,i1,i2)
               end do
            end do
         end do
        end do

        do i=1,nparts      ! symmetrize h along diagonal
          do l=2,ndim
          do m=1,l-1
             h1(m,l,i,i)=h1(l,m,i,i)
             h2(m,l,i,i)=h2(l,m,i,i)
          enddo
          enddo
        enddo


        xi=0.d0
        xa=0.d0
        do i=1,nparts
         it=idtype(i)
         ehb=hbs2m(it)
         do l=1,ndim
            xi=xi+force1(l,i)*force2(l,i)
            xa=xa+force1(l,i)*dforce2(l,i)
     .           +dforce1(l,i)*force2(l,i)
            xdummy=force1(l,i)*h2(l,l,i,i)
     .            +force2(l,i)*h1(l,l,i,i)
            xa=xa+2.d0*h1(l,l,i,i)*h2(l,l,i,i)*ehb
            do m=1,l-1
              xdummy=xdummy
     .        +force1(m,i)*h2(l,m,i,i)
     .        +force2(m,i)*h1(l,m,i,i)
              xa=xa+4.d0*h1(l,m,i,i)*h2(l,m,i,i)*ehb
            end do
            do m=l+1,ndim
              xdummy=xdummy
     .        +force1(m,i)*h2(l,m,i,i)
     .        +force2(m,i)*h1(l,m,i,i)
            end do
            gradu(l,i)=gradu(l,i)+plambdat*xdummy
         end do
       end do
       u=u+plambdat*xi
       del2u=del2u+plambdat*xa
       xa=0.d0
       xi=0.d0
       xa2=0.d0
       do i1=1,nparts
        do i2=1,i1-1
          it=idtype(i1)
          jt=idtype(i2)
          ehb=hbs2m(it)+hbs2m(jt)
          do l=1,ndim
            xd1=0.d0
            xd2=0.d0
            xi=xi+v(l,i1,i2)*w(l,i1,i2) ! 2-body due to 3-body
            xa2=xa2+d2v(l,i1,i2)*w(l,i1,i2)
     .             +d2w(l,i1,i2)*v(l,i1,i2) ! 2-body due to 3-body
            do m=1,ndim
              xd1=xd1+force1(m,i2)*h2(m,l,i1,i2)
     .               +force2(m,i2)*h1(m,l,i1,i2)
              xd1=xd1+1.d0*(h1(l,m,i1,i2)*w(m,i1,i2)
     .                     +h2(l,m,i1,i2)*v(m,i1,i2)) ! 2-body due to 3-body
              xd2=xd2+force1(m,i1)*h2(m,l,i1,i2)
     .               +force2(m,i1)*h1(m,l,i1,i2)
              xd2=xd2-1.d0*(h1(l,m,i1,i2)*w(m,i1,i2)
     .                     +h2(l,m,i1,i2)*v(m,i1,i2)) ! 2-body due to 3-body
              xa=xa+2.d0*ehb
     .           *h1(m,l,i1,i2)*h2(m,l,i1,i2)
              xa2=xa2+2.d0*h1(m,l,i1,i2)*h2(m,l,i1,i2)*ehb ! 2-body due to 3-body
            end do
            gradu(l,i1)=gradu(l,i1)+plambdat*xd1
            gradu(l,i2)=gradu(l,i2)+plambdat*xd2
          end do
        end do
       end do
       u=u-1.d0*plambdat*xi  ! 2-body due to 3-body
       del2u=del2u+plambdat*xa
     .          -1.d0*plambdat*xa2
      end if



! mh end


      return
      end