******************************************************************************
C                          SUBROUTINE PERM.F
******************************************************************************
C This version works best for argcut<8.0.  For larger values, the original
C program runs faster.  The speed of the original version
C is insensitive to the value of argcut.  
******************************************************************************
      subroutine perm(x,y,z,t,matoms,iatom,iatomp,sum,con,
     x                islice,iptot)
      include 'pimc_sqr3sc.dcl'
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension t(nperm),f(nat,nat),sum(0:nperm)
      integer ikeep(nfact+nat),jkeep(nfact+nat)
      integer matoms(nperm),iatom(nperm,4),iatomp(nperm,4)

*//////////////////////////////////////////////////////////////////////////////
C                         ONE- AND TWO-BODY FACTORS
C Here we calculate the two-particle factors f(ip), whose products will be
C used to construct two-,three-, and four-particle permutations.  See equation
C (49).  There are np2 values of f(ip), where np2=nat*(nat-1).  
C
C Do not calculate these factors for static atoms.  Note also that static
C atoms are always distinguishable, so its only necessary to check the static
C flag at the first time slice.  
*//////////////////////////////////////////////////////////////////////////////

      icount1=0
      icount2=0
      icount3=0
      icount4=0
      sum1=0.
      sum2=0.
      sum3=0.
      sum4=0.

      cutoff=exp(-argcut)
      ip=0
      isl2=islice+nbis            ! Index all 1-body perms first
      do i=1,nat
         f(i,i)=0.0
         call dstpbc(x(islice,i),x(isl2,i),xdim,xd1,sgn)
         call dstpbc(y(islice,i),y(isl2,i),ydim,yd1,sgn)
         zd1=z(islice,i)-z(isl2,i)
         arg=con*(xd1**2+yd1**2+zd1**2)
         ip=ip+1
         ikeep(ip)=i
         jkeep(ip)=i
         f(i,i)=exp(-arg)
      end do
      np1=ip
      do i=1,nat
         do j=1,nat
            if(i.eq.j) goto 1
            f(i,j)=0.0
            call dstpbc(x(islice,i),x(isl2,j),xdim,xd2,sgn)
            call dstpbc(y(islice,i),y(isl2,j),ydim,yd2,sgn)
            zd2=z(islice,i)-z(isl2,j)
            arg=con*(xd2**2+yd2**2+zd2**2)
            if(arg.lt.argcut) then
               ip=ip+1
               ikeep(ip)=i
               jkeep(ip)=j
               f(i,j)=exp(-arg)
            end if
 1          continue
         end do
      end do
      np2=ip-np1
      
*////////////////////////////////////////////////////////////////////////////
C                          IDENTITY PERMUTATION
*////////////////////////////////////////////////////////////////////////////

      iperm=0
      sum(iperm)=0.0
      do i=1,np1
         icount1=icount1+1
         ii=ikeep(i)
         iperm=iperm+1
         t(iperm)=f(ii,ii)
         sum1=sum1+t(iperm)
         sum(iperm)=sum(iperm-1)+t(iperm)
         matoms(iperm)=1
         iatom(iperm,1)=ii
         iatomp(iperm,1)=ii
      end do

*//////////////////////////////////////////////////////////////////////////////
C                        MULTI-BODY PERMUTATIONS
C Two-, three-, and four-body cyclic permutations.  Not all 4-body perms 
C are given.
*//////////////////////////////////////////////////////////////////////////////

      do i=np1+1,np2
         ii=ikeep(i)
         jj=jkeep(i)
                                ! 2-body, ij-->ji
         prod=f(ii,jj)*f(jj,ii)
         if(prod.gt.cutoff)then
            icount2=icount2+1
            iperm=iperm+1
            t(iperm)=prod
            sum2=sum2+t(iperm)
            sum(iperm)=sum(iperm-1)+t(iperm)
            matoms(iperm)=2
            iatom(iperm,1)=ii
            iatom(iperm,2)=jj
            iatomp(iperm,1)=jj
            iatomp(iperm,2)=ii
         end if
         do j=np1+1,np2
            if(j.eq.i) goto 10  ! So ii.ne.kk and jj.ne.ll
            kk=ikeep(j)
            ll=jkeep(j)
                                ! 3-body, ijk-->kij
            if(ll.eq.ii.and.kk.ne.jj) then
               prod=f(ii,jj)*f(jj,kk)*f(kk,ll)
               if(prod.gt.cutoff)then
                  icount3=icount3+1
                  iperm=iperm+1
                  t(iperm)=prod
                  sum3=sum3+t(iperm)
                  sum(iperm)=sum(iperm-1)+t(iperm)
                  matoms(iperm)=3
                  iatom(iperm,1)=ii
                  iatom(iperm,2)=jj
                  iatom(iperm,3)=kk
                  iatomp(iperm,1)=kk
                  iatomp(iperm,2)=ii
                  iatomp(iperm,3)=jj
               end if
            end if
                                 ! 3-body, ijk-->jki
            if(ll.eq.jj.and.kk.ne.ii) then
               prod=f(ii,kk)*f(jj,ii)*f(kk,jj)
               if(prod.gt.cutoff) then
                  icount3=icount3+1
                  iperm=iperm+1
                  t(iperm)=prod
                  sum3=sum3+t(iperm)
                  sum(iperm)=sum(iperm-1)+t(iperm)
                  matoms(iperm)=3
                  iatom(iperm,1)=ii
                  iatom(iperm,2)=jj
                  iatom(iperm,3)=kk
                  iatomp(iperm,1)=jj
                  iatomp(iperm,2)=kk
                  iatomp(iperm,3)=ii
               end if
            end if
                                ! 4-body, ijkl-->lijk  
            if(jj.ne.kk.and.ii.ne.ll.and.jj.ne.ll.and.ii.ne.kk)then
               prod=f(ii,jj)*f(jj,kk)*f(kk,ll)*f(ll,ii)
               if(prod.gt.cutoff) then
                  icount4=icount4+1
                  iperm=iperm+1
                  t(iperm)=prod
                  sum4=sum4+t(iperm)
                  sum(iperm)=sum(iperm-1)+t(iperm)
                  matoms(iperm)=4
                  iatom(iperm,1)=ii
                  iatom(iperm,2)=jj
                  iatom(iperm,3)=kk
                  iatom(iperm,4)=ll
                  iatomp(iperm,1)=ll
                  iatomp(iperm,2)=ii
                  iatomp(iperm,3)=jj
                  iatomp(iperm,4)=kk
               end if
                                ! 4-body, ijkl-->jkli
               prod=f(ii,ll)*f(jj,ii)*f(kk,jj)*f(ll,kk)
               if(prod.gt.cutoff) then
                  icount4=icount4+1
                  iperm=iperm+1
                  t(iperm)=prod
                  sum4=sum4+t(iperm)
                  sum(iperm)=sum(iperm-1)+t(iperm)
                  matoms(iperm)=4
                  iatom(iperm,1)=ii
                  iatom(iperm,2)=jj
                  iatom(iperm,3)=kk
                  iatom(iperm,4)=ll
                  iatomp(iperm,1)=jj
                  iatomp(iperm,2)=kk
                  iatomp(iperm,3)=ll
                  iatomp(iperm,4)=ii
               end if
                                ! 4-body, ijkl-->klij
               prod=f(ii,kk)*f(jj,ll)*f(kk,ii)*f(ll,jj)
               if(prod.gt.cutoff) then
                  icount4=icount4+1
                  iperm=iperm+1
                  t(iperm)=prod
                  sum4=sum4+t(iperm)
                  sum(iperm)=sum(iperm-1)+t(iperm)
                  matoms(iperm)=4
                  iatom(iperm,1)=ii
                  iatom(iperm,2)=jj
                  iatom(iperm,3)=kk
                  iatom(iperm,4)=ll
                  iatomp(iperm,1)=kk
                  iatomp(iperm,2)=ll
                  iatomp(iperm,3)=ii
                  iatomp(iperm,4)=jj
               end if
            end if
 10         continue
         end do
      end do
      iptot=iperm
c      print*,sum1/sum(iperm),sum2/sum(iperm),sum3/sum(iperm),
c     x     sum4/sum(iperm)
c      print*,icount1,icount2,icount3,icount4
c      print*
      return
      end 
