program hard_spheres_mc
integer i,j,k,t,tmax,p,n;  parameter(p=5,n=p**3,tmax=100*n)
real r(1:3,1:n),rt(1:3),rij2(1:n),u(1:4),d,delta; parameter(d=0.5,delta=0.1)

forall(k=1:3,i=1:n) r(k,i)=modulo((i-1)/p**(k-1),p)+0.5

do t=1,tmax
   call random_number(u);   i=int(n*u(1))+1
   rt=r(:,i)+delta*(2*u(2:4)-1)
   if ( (minval(rt)>d/2).and.(maxval(rt)<p-d/2)) then
      forall(j=1:n,j/=i) rij2(j)=sum((rt-r(:,j))**2); rij2(i)=2*d**2
      if (minval(rij2)>d**2) then
         r(:,i)=rt ; write(6+i,*) r(:,i)
         endif; endif
enddo
end

! gnuplot>  splot for [i=7:131]  'fort.'.i with lines notitle
