***********************************************************************
C                           TRIANGLE_GRSOL.F
C Initializes helium positions on a triangular lattice above the 
C substate.
***********************************************************************
      subroutine initialize(x,y,z,xp,yp,zp)
      include 'pimc_sqr3sc.dcl'
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension xp(2*m,nat),yp(2*m,nat),zp(2*m,nat)

      idum=-78915
      inc=0
      
      den=nat/xdim/ydim
c      a1=sqrt(sqrt(3.)/den/2.)
      a1=4.26
      a2=sqrt(3.)*a1
      zstart=3.2
      zstep=sig

      rx1=a1/2.
      ry1=a2/2.

      del=.2*sig
      xstart=del+2.46/2./sqrt(3.)
      ystart=del+1.23   

      incx=1
      incy=1
      incz=1
      
      do j=1,2*(nat/2),2
         do islice=1,m
            x(islice,j)=xstart+a1*(incx-1)+del*(.5-ran1(idum))
            y(islice,j)=ystart+a2*(incy-1)+del*(.5-ran1(idum))
            z(islice,j)=zstart+zstep*(incz-1)+del*(.5-ran1(idum))
            x(islice,j+1)=xstart+a1*(incx-1)+rx1+del*(.5-ran1(idum))
            y(islice,j+1)=ystart+a2*(incy-1)+ry1+del*(.5-ran1(idum))
            z(islice,j+1)=zstart+zstep*(incz-1)+del*(.5-ran1(idum))

            xp(islice,j)=x(islice,j)
            yp(islice,j)=y(islice,j)
            zp(islice,j)=z(islice,j)
            x(islice+m,j)=x(islice,j)
            y(islice+m,j)=y(islice,j)
            z(islice+m,j)=z(islice,j)
            xp(islice+m,j)=xp(islice,j)
            yp(islice+m,j)=yp(islice,j)
            zp(islice+m,j)=zp(islice,j)

            xp(islice,j+1)=x(islice,j+1)
            yp(islice,j+1)=y(islice,j+1)
            zp(islice,j+1)=z(islice,j+1)
            x(islice+m,j+1)=x(islice,j+1)
            y(islice+m,j+1)=y(islice,j+1)
            z(islice+m,j+1)=z(islice,j+1)
            xp(islice+m,j+1)=xp(islice,j+1)
            yp(islice+m,j+1)=yp(islice,j+1)
            zp(islice+m,j+1)=zp(islice,j+1)
         end do
         incx=incx+1
         if(xstart+a1*(incx-1)+rx1.gt.xdim) then
            incx=1
            incy=incy+1
            if(ystart+a2*(incy-1)+ry1.gt.ydim) then
               incy=1
               incz=incz+1
            end if
         end if
      end do
      if(2*(nat/2).ne.nat) then
         xran=xdim*ran1(idum)
         yran=ydim*ran1(idum)
         do islice=1,m
            x(islice,nat)=xstart+a1*(incx-1)+del*(.5-ran1(idum))
            y(islice,nat)=ystart+a2*(incy-1)+del*(.5-ran1(idum))
            z(islice,nat)=zstart+zstep*(incz-1)+del*(.5-ran1(idum))

            xp(islice,nat)=x(islice,nat)
            yp(islice,nat)=y(islice,nat)
            zp(islice,nat)=z(islice,nat)
            x(islice+m,nat)=x(islice,nat)
            y(islice+m,nat)=y(islice,nat)
            z(islice+m,nat)=z(islice,nat)
            xp(islice+m,nat)=xp(islice,nat)
            yp(islice+m,nat)=yp(islice,nat)
            zp(islice+m,nat)=zp(islice,nat)
         end do
      end if
      return
      end
      
