implicit double precision(a-h,o-z)
      real*8 x(3,110000),vx(330000),vy(110000),vz(110000)
      data wkb/1.380658d-23/
      open(3,file='xe_ini_1.d',status='unknown')
c input
      sigma = 4.100d-10
      e     = 0.305d-20
      wm    = 2.1802d-25
      nlat  = 20
      npart = 4*nlat**3
      rho   = 750.d0

      const = 6.0186d-9
      cell  = const*(npart/rho)**(1.d0/3.d0)

      tset  = 300.d0
      dt    = 0.5d-14
      nstep = 10000

      wnp     = (npart/4.0d0)**(1.0d0/3.0d0) - 1.0d0
      a       = cell / (wnp+1.0d0)
      np=nint(wnp)
      ijk = 0
      do lg=0,1
        do i=0,np
          do j=0,np
            do k=0,np
              ijk =ijk + 1
              x(1,ijk) = i * a + lg * a * 0.5d0
              x(2,ijk) = j * a + lg * a * 0.5d0
              x(3,ijk) = k * a
            enddo
          enddo
        enddo
      enddo
      do lg=1,2
        do i=0,np
          do j=0,np
            do k=0,np
              ijk = ijk + 1
              x(1,ijk) = i * a + (2-lg) * a * 0.5d0
              x(2,ijk) = j * a + (lg-1) * a * 0.5d0
              x(3,ijk) = k * a + a * 0.5d0
            enddo
          enddo
        enddo
      enddo
      do i=1,npart
        vx(npart+i)=vy(i)
        vx(2*npart+i)=vz(i)
      enddo
      call nrml2(tset,vx,3*npart,wm,wkb)
      do i=1,npart
        vy(i)=vx(npart+i)
        vz(i)=vx(2*npart+i)
      enddo
c ---- output of the initial conditions ----
      time=0.d0
      write(3,3456) npart,nstep,sigma,e,wm,cell,tset,dt
 3456 format(2i6,6d23.15)
      do k=1,npart
        write(3,2001) x(1,k),x(2,k),x(3,k)
        write(3,2001) vx(k),vy(k),vz(k)
      enddo
      write(3,3004) time
 2001 format(3d23.15)
 3001 format(2d23.15)
 3002 format(3d23.15,d23.15,i8)
 3004 format(d23.15)
 4000 format(3d23.15/2d23.15,i8)
      stop
      end
c======================================================================
      subroutine nrml2(tset,a,n,wm,wkb)
      implicit double precision(a-h,o-z)
      real*8 a(330000),wm,wkb

      data b,c,d,e/1.256853d0,1.706340d0,-0.2898778d0,-0.05971928d0/
      data p,q,r/1.451153d0,-0.3040790d0,-0.5111851d0/

      sd=dsqrt(wkb*tset/wm)
      am=0.0d0
      iseed=1314
      call unfm2(iseed,a,n)
      do i=1,n
        u=a(i)
        x=u-0.5d0
        t=x*x
        if(dabs(x).le.0.46875d0) then
          g=(e/(t+d)+t+b)*c*x
        else
          v=sqrt(-dlog(0.5d0-dabs(x)))
          g=dsign((r/v+q+v)*p,x)
        endif
        a(i)=g*sd+am
      enddo
      return
      end
c======================================================================
      subroutine unfm2(ix,vh,n)
      implicit double precision(a-h,o-z)
      real*8 vh(330000)
      real*8 aa,b,c,ci,x
      data aa/32771.d0/,b/1234567891.d0/
      data c/2147483648.d0/,ci/4.6566128730773925d-10/
      x=ix
      do i=1,n
        x = dmod( aa * x + b,c)
        vh(i) = x * ci
      enddo
      ix = x
      return
      end