c****************************************************************************
c	This is the source code for simulating symmetric non-additive 
c	hard spheres. The program takes as input the number of molecules,
c	simulation time and non additivity constant.
c	The initial conditions are given as the de-mixed closed pack 
c	hexagonal packing.
c	The density is decreased by a factor after nstep moves.
c	This is repeated by the loop "jjj".
c	One has to stop the program when desired denstiy range is simulated.
c	The program uses a neighbor list. The neighbor list is updated after 
c	each move.
c	Each particle can either swap or move at each step.
c	n:number of spheres
c	nstep: simulation time
c	fac: non additivity factor
c*****************************************************************************
	program main
	real*8 den,fac
	integer nstep,n
	open(131,file='inpu')
	nstep=10000
	n=44**2
	
	read(131,*) fac

	write(*,*) fac
	call hards(n,nstep,fac)
	end


	subroutine hards(n,nstep,fac)
	integer n,h,box(9),sak1,sak2,sak3,sak4,swaprej
	real*8 lx,ly,den,rann,fac
	parameter(mmax=400)
	integer type(n),rej,oldicell
	integer list(0:n+mmax**2),nbox(9,mmax**2)
	real*8 x(n,2),ran(2),xn(2),sg(0:1),sig,del
	open(10,file='out')
	open(11,file='celllen')
	write(11,*) 'fac,n,nstep,20===',fac,n,nstep
	nbox=0
	list=0
c*****************************************************************************
c	sg: diameters of the spheres
c	cof: length scaling factor for decreasing density
c	del: initial Metropolis step size. With in the program this step size is 
c	adjusted to satisfy acceptance ratio between 0.4-0.6
c*****************************************************************************
	cof=1.001
	sg(1)=2.d0
	sg(0)=2.d0
	del=sg(1)/2.d0

	lx=dsqrt(dble(n))*sg(1)
	ly=dsqrt(0.75d0*n)*sg(1)+2.2d0*(fac-1.d0)*dsqrt(0.75d0)*sg(1)
	den=n*(sg(1)**2.d0)/lx/ly
	write(11,*) lx,ly,den,den*fac**2.d0
	k=ceiling(dsqrt(dble(n)))

	mx=floor(lx/(sg(1)*fac))
	my=floor(ly/(sg(1)*fac))
	celx=lx/mx
	cely=ly/my
	lx=celx*dble(mx)
	ly=cely*dble(my)
	write(*,*) 'mx',mx,my
	j=0
	jj=0
c*****************************************************************************
c	initial configuration
c*****************************************************************************
	do i=1,n
		x(i,1)=lx/k*j+0.0001
		if(mod(jj,2).eq.1) x(i,1)=lx/k*(j+.5d0)
		x(i,2)=sg(1)*jj*dsqrt(0.75d0)+0.0001
		if(j.lt.k-1) then
			j=j+1
		else
			j=0
			jj=jj+1
		endif
	enddo
	x(n/2+1:n,2)=x(n/2+1:n,2)+(ly-x(n,2)-sg(1)*dsqrt(.75d0))/2.d0
	do i=1,n/2
	type(i)=0
	type(i+n/2)=1
	enddo

	do i=1,n
		write(10,20) type(i),x(i,:)
	enddo
	
20	format(i3,2f10.5)
	write(10,*)

c*****************************************************************************
c	Create  neighbor list
c*****************************************************************************
	do jj=1,n
		icell=ceiling(x(jj,1)/celx)+(ceiling(x(jj,2)/cely)-1)*mx
		list(jj)=list(n+icell)
		list(n+icell)=jj
	enddo

	do i=1,mx*my
		call check(i,mx,my,nbox(:,i))
	enddo

	do jjj=1,1000
	iswap=0
	imove=0
	swaprej=0
	rej=0
	x=x*cof
	celx=celx*cof
	cely=cely*cof
	lx=celx*dble(mx)
	ly=cely*dble(my)
	den=den/(cof**2.d0)
	write(11,*) lx,ly,den,den*fac**2.d0
	do 1000 ii=1,nstep
		do i=1,n
			rann=0.d0
			call random_number(rann)
c********************************************************
c	SWAP
c********************************************************
			if(rann.ge.0.5d0) then
				iswap=iswap+1

				call random_number(rann)
				irann=ceiling(rann*n/2)
				if(i.le.n/2) then
					irann=irann+n/2
				endif
				icell=ceiling(x(i,1)/celx)+(ceiling(x(i,2)/cely)-1)*mx
				iicell=icell
			j=0
130			j=j+1
			if(j.gt.9) then
			goto 233
			endif
			nbr=list(n+nbox(j,icell))
131			if(nbr.eq.i) nbr=list(nbr)
			if(nbr.eq.0) goto 130
			disx=abs(x(i,1)-x(nbr,1))
			if(disx.gt.(lx/2.)) disx=lx-disx
			disy=abs(x(i,2)-x(nbr,2))
			if(disy.gt.(ly/2.d0)) disy=ly-disy
			dis=dsqrt(disx**2.d0+disy**2.d0)
			sig=(sg(type(irann))+sg(type(nbr)))/2.d0
			if(type(irann).ne.type(nbr)) sig=sig*fac
			if(dis.gt.sig) then
				nbr=list(nbr)
				goto 131
			else
				swaprej=swaprej+1
				goto 40
			endif

233			icell=ceiling(x(irann,1)/celx)+(ceiling(x(irann,2)/cely)-1)*mx
			j=0
230			j=j+1
			if(j.gt.9) then
			if(iicell.eq.icell) then
			xn=x(i,:)
			x(i,:)=x(irann,:)
			x(irann,:)=xn
			goto 40
			endif

			kkk=n+nbox(1,icell)
234			if(list(kkk).eq.irann) then
			sak1=kkk
			sak2=list(list(kkk))
			else
			kkk=list(kkk)
			goto 234
			endif
				icell=iicell
			kkk=n+nbox(1,icell)
235			if(list(kkk).eq.i) then
			sak3=kkk
			sak4=list(list(kkk))
			else
			kkk=list(kkk)
			goto 235
			endif
			
			xn=x(i,:)
			x(i,:)=x(irann,:)
			x(irann,:)=xn
			list(sak1)=i
			list(i)=sak2
			list(sak3)=irann
			list(irann)=sak4			

			goto 40
			endif


			nbr=list(n+nbox(j,icell))
231			if(nbr.eq.irann) nbr=list(nbr)
			if(nbr.eq.0) goto 230
			
			disx=abs(x(irann,1)-x(nbr,1))
			if(disx.gt.(lx/2.d0)) disx=lx-disx
			disy=abs(x(irann,2)-x(nbr,2))
			if(disy.gt.(ly/2.d0)) disy=ly-disy
			dis=dsqrt(disx**2.d0+disy**2.d0)
			sig=(sg(type(i))+sg(type(nbr)))/2.d0
			if(type(i).ne.type(nbr)) sig=sig*fac
	
			if(dis.gt.sig) then
				nbr=list(nbr)
				goto 231
			else
				swaprej=swaprej+1
				goto 40
			endif



	
c********************************************************
c	MOVE
c********************************************************
			else
		imove=imove+1
3001			call random_number(ran)
			ran=ran*del-del/2.d0
			xn=x(i,:)-ran	
			if(xn(1).ge.lx) then
				xn(1)=xn(1)-lx
			elseif(xn(1).lt.0.d0) then
				xn(1)=lx+xn(1)
			endif			
			if(xn(2).ge.ly) then
				xn(2)=xn(2)-ly
			elseif(xn(2).lt.0.d0) then
				xn(2)=ly+xn(2)
			endif
			icell=ceiling(xn(1)/celx)+(ceiling(xn(2)/cely)-1)*mx

			j=0
30			j=j+1
			if(j.gt.9) then
				oldicell=ceiling(x(i,1)/celx)+(ceiling(x(i,2)/cely)-1)*mx	
				x(i,:)=xn(:)
				if(oldicell.eq.icell) goto 40

				kk=list(n+oldicell)
				kkk=n+oldicell
32				if(kk.eq.i) then
					list(kkk)=list(i)
				else
					kkk=kk
					kk=list(kkk)
					goto 32
				endif
				
				list(i)=list(n+icell)
				list(n+icell)=i


				goto 40
			endif
			nbr=list(n+nbox(j,icell))
31			if(nbr.eq.i) nbr=list(nbr)
			if(nbr.eq.0) then
				goto 30
			endif
			disx=abs(xn(1)-x(nbr,1))
			if(disx.gt.(lx/2.d0)) disx=lx-disx
			disy=abs(xn(2)-x(nbr,2))
			if(disy.gt.(ly/2.d0)) disy=ly-disy
			dis=dsqrt(disx**2.d0+disy**2.d0)
			sig=(sg(type(i))+sg(type(nbr)))/2.d0
			if(type(i).ne.type(nbr)) sig=sig*fac
	
			if(dis.gt.sig) then
				nbr=list(nbr)
				goto 31
			else
				rej=rej+1
			endif
	endif
40		enddo	
		if(float(rej)/imove.lt.0.4d0) then
		del=del*1.1
		write(*,*) 'increased',del
		elseif(float(rej)/imove.gt.0.6d0) then 
		if(del.gt.0.1) then
		del=del/1.1
		write(*,*) 'decreased',del
		endif
		endif

		if(mod(ii,1000).eq.0) write(*,*) ii,float(rej)/imove,float(swaprej)/iswap
		if(mod(ii,50).eq.0) then
			do i=1,n
				write(jjj+9000,20) type(i),x(i,:)
			enddo
			write(jjj+9000,20) 
		endif
1000	enddo
	enddo
	write(*,*) rej,swaprej,float(rej)/imove,float(swaprej)/iswap
	return
	end


	Subroutine check(icell,mx,my,box)
c*****************************************************************************
c	For the neighbor list create the box mesh
c*****************************************************************************

	integer box(9),row,col,icell,mx,my,rowm,rowp,colm,colp
	row=ceiling(float(icell)/float(mx))
	col=icell-(row-1)*mx
	if(row.eq.1) then
	rowm=my
	rowp=2
	elseif(row.eq.my) then
	rowm=my-1
	rowp=1
	else
	rowm=row-1
	rowp=row+1
	endif
	if(col.eq.1) then
	colm=mx
	colp=2
	elseif(col.eq.mx) then
	colm=mx-1
	colp=1
	else
	colm=col-1
	colp=col+1
	endif
	box(1)=icell
	box(2)=mx*(rowm-1)+col
	box(3)=mx*(rowp-1)+col
	box(4)=mx*(row-1)+colm
	box(5)=mx*(row-1)+colp
	box(6)=mx*(rowm-1)+colm
	box(7)=mx*(rowm-1)+colp
	box(8)=mx*(rowp-1)+colm
	box(9)=mx*(rowp-1)+colp

	return
	end

