c COMPILE USING F77, e.g.  f77 -o xorder ordering.f
c EXECUTABLE is then XORDER, and may be ran by typing  ./xorder
c By type ./xorder > file.name, you may save unit 5 (i.e. screen) output.
c
c
	program main
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	Chemical ordering in a two-dim square lattice
c			by diffusion of one vacancy
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	implicit none
	integer*4 ns,noutmx,nssq
	parameter (ns=30,noutmx=1000)   ! if ns > 50, change last write format
	parameter (nssq=ns*ns)
	integer*4 site(ns,ns),nn(ns,ns),nout,ninit
	integer*4 modl(-1:ns+2),vacx,vacy,nvec(2,4)
	integer*4 i,ii,j,jj,k,kk,x(4),y(4),natot,nbtot
	integer*4 iseed,nstep
	real*8 vaa,vab,vbb,temp,tempk,bconc,RP_BA
	real*8 texpaa(0:3),texpab(0:3),texpbb(0:3)
	real*8 sro(noutmx),lro(noutmx),time(noutmx),atime
	real*8 w(4),dtime,wcumul,xx,yy
	real*8 ranno,ran
	data site,nn /nssq*1, nssq*0/
	data sro,lro,time /noutmx*0. , noutmx*0., noutmx*0./
	iseed=1234567
c
ccccccccccccccccccccccc  VARIABLE   DETAILS 
c	ns : number of lattice sites in x or y direction
c	nssq = ns*ns : total number of lattice sites
c	nstep : total number of vacancy jump
c	nout : number of sro and lro measurements
c	noutmx : maximum number of nout
c
c	ninit : if equal 1, build an ordered initial configuration;
c		otherwise build it randomly
c
c	tempk : temperature in Kelvin
c
c	bconc : concentration of B species
c
c	vaa,vab,vbb : interaction energies for AA,AB,BB pairs
c
c	site(ns,ns) : this array gives the occupation of any site
c			1 for an A atom
c			0 for the vacancy
c	      	       -1 for a B atom
c	vacx,vacy : x- and y-coordinate of the vacancy
c
c	nn(ns,ns) : this array gives the number of A nearest neighbors
c			of the considered site
c
c	modl(-1:ns+2) : returns the input value modulo ns, 
c			except modl(ns)=ns
c
c	nvec(2,4) : gives the x and y components of the four nearest
c			neighbor vectors
c
c	texpaa,texpab,texpbb : tabulated exponentials
c
c	w(4) : array of jump frequencies
c	wcumul : sum of the four w(i)
c
c	dtime : time spent between two jumps
c	atime : time elapsed from the very first jump
c
c	time(nout) : current value of the total time
c
c	sro(nout),lro(nout) : current values of Short Range Order
c		and Long Range Order parameters
c
c	Output : code produces four files  
c           ordera.res  w/ locations of A-atoms
c           orderb.res  w/ locations of B-atoms
c           orderv.res  w/ locations of vacancy
c            order.res  w/ 2-dim listing of 1,0,-1  (A, Vac, B, resp)
c
c	Could use Kaleidagraph to draw  this final atomic configuration
c       or just look at order.res
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccc
c	Reading data
ccccccccccccccccccccccc
	open (unit=5,file='order.data',status='old')
	read (5,*) nstep,nout,ninit
	read (5,*) tempk,bconc
	read (5,*) vaa,vab,vbb
	temp=tempk/11604.9
	close (5)
ccccccccccccccccccccccc
c	Initialization
ccccccccccccccccccccccc
c	tabulated exponentials
	do 10 i=0,3
	  texpaa(i)=exp(i*vaa/temp)
	  texpab(i)=exp(i*vab/temp)
	  texpbb(i)=exp(i*vbb/temp)
 10	continue
c--------------------------------
c	Modulo ns table
c--------------------------------
	do 20 i=-1,ns+2
	  modl(i)=mod(i,ns)
	  if (modl(i).le.0) modl(i)=modl(i)+ns
 20	continue
c--------------------------------
c	Array of nearest neighbor vectors
c-------------------------------
	nvec(1,1)=1
	nvec(2,1)=0
	nvec(1,2)=-1
	nvec(2,2)=0
	nvec(1,3)=0
	nvec(2,3)=1
	nvec(1,4)=0
	nvec(2,4)=-1
c---------------------------------
c	Building the initial configuration
c----------------------------------
	nbtot=nint(bconc*real(ns)**2)
	natot=ns**2-nbtot
	if (ninit.ne.1) then 
	  do 30 i=1,nbtot
 35	    continue
	    j=int(ranno(iseed)*real(ns))+1
	    jj=int(ranno(iseed)*real(ns))+1
	    if (site(j,jj).eq.1) then
	      site(j,jj)=-1
	    else
	      goto 35
	    end if
 30	  continue
	else
	  k=0
	  do 36 j=1,ns
	    do 37 i=1,ns,2
	      if (mod(j,2).eq.0) then
	        site(i+1,j)=-1
	      else
	        site(i,j)=-1
	      end if
	    k=k+1
	    if (k.eq.nbtot) goto 38
 37	    continue
 36	  continue
	end if
 38	continue
	j=int(ranno(iseed)*real(ns))+1
	jj=int(ranno(iseed)*real(ns))+1
	site(j,jj)=0
	vacx=j
	vacy=jj
c----------------------------------------
c	Computing the number of A nearest neighbors of any site
c----------------------------------------
	do 40 i=1,ns
	  do 50 j=1,ns
	    do 60 k=1,4
	    ii=modl(i+nvec(1,k))
	    jj=modl(j+nvec(2,k))
	    if (site(ii,jj).eq.1) nn(i,j)=nn(i,j)+1
 60	    continue
 50	  continue
 40	continue
c---------------------------------------
c	Calculating initial sro and lro parameters
c---------------------------------------
c
c------------------------------------------------------------------------------
c DEFINITIONS OF SRO CAN BE ARBITRARY, sometimes just for convienience.
c
c STANDARD: The usual SRO (i.e., Warren-Cowley) = ALPHA_ij =   1 -  P_BA/P_A
c  where P_BA = joint prob. of an A as neighbor given a B at central site i.
c  and   P_A  = prob. of A at a site j, which is composition of A.
c  NOTE: with this definition sro-random= 0, sro-phase-sep=+1, and sro-ord= -1
c        The LRO is always between  -1 < LRO < 1
c
c POSSIBLE: SRO= P_B*P_BA - Z*P_B*(1-P_B), where Z=4 on square lattice.
c  Among other things this mean that you do not blow-up if you input bconc=0.
c  NOTE: with this definition sro-random= 0, sro-phase-sep=-1, and sro-ord= +1
c        The LRO is always between  -1 < LRO < 1
c        SRO then mirrors LRO but it is not direclty related to probability
c------------------------------------------------------------------------------
c
	  time(1)=0.
	  do 240 ii=1,ns
	    do 250 jj=1,ns
	      if (site(ii,jj).lt.0) sro(1)=sro(1)+real(nn(ii,jj))
 250	    continue
 240	  continue
c
cProbability of B-A neighbors for Random Alloy w/ B atom central site
          RP_BA= 4.*bconc*(1.-bconc)
c
          sro(1)=sro(1)/(real(ns)**2)            ! P_B*P_BA
cNon-standard
c         sro(1)=sro(1)- RP_BA                   ! different normalization, c=0 OK
cWCdef    
          sro(1)= 1 - sro(1)/RP_BA               ! 1 - P_BA/P_A 
	  do 260 jj=1,ns
	    do 270 ii=1,ns,2
	      if (mod(jj,2).eq.0) then
	        if (site(ii+1,jj).eq.1) lro(1)=lro(1) +1.
	      else
	        if (site(ii  ,jj).eq.1) lro(1)=lro(1) +1.
	      end if
 270	    continue
 260	  continue
	  lro(1)=abs(2.*lro(1)/real(natot)-1.)
ccccccccccccccccccccccc
c	Monte-Carlo loop
ccccccccccccccccccccccc
	atime=0.
	dtime=0.
	nstep=nstep/nout
	do 100 i=1,nout
	  do 110 j=1,nstep
c---------------------------------------
c	Computing the four jump frequencies
c---------------------------------------
	  do 115 k=1,4
	  x(k)=modl(vacx+nvec(1,k))
	  y(k)=modl(vacy+nvec(2,k))
	  if (site(x(k),y(k)).eq.1) then
	    w(k)=texpaa(nn(x(k),y(k)))*texpab((3-nn(x(k),y(k))))
	  else
	    w(k)=texpab(nn(x(k),y(k)))*texpbb((3-nn(x(k),y(k))))
	  end if
 115	  continue
c--------------------------------------
c	Choosing the jump to be performed
c--------------------------------------
	  dtime=1/(w(1)+w(2)+w(3)+w(4))
	  wcumul=0.
	  ran=ranno(iseed)/dtime
	  do 120 k=1,4
	    wcumul=wcumul+w(k)
	    if (ran.le.wcumul) goto 125
 120	  continue
 125	  continue
	  atime=atime+dtime
c--------------------------------------
c	Updating neighborhoods and positions
c--------------------------------------
	  if (site(x(k),y(k)).eq.1) then
	    do 130 kk=1,4
	      xx=modl(x(k)+nvec(1,kk))
	      yy=modl(y(k)+nvec(2,kk))
	      nn(xx,yy)=nn(xx,yy)-1
 130	    continue
	    do 135 kk=1,4
	      xx=modl(vacx+nvec(1,kk))
	      yy=modl(vacy+nvec(2,kk))
	      nn(xx,yy)=nn(xx,yy)+1
 135	    continue
	  end if
	  site(vacx,vacy)=site(x(k),y(k))
	  vacx=x(k)
	  vacy=y(k)
	  site(vacx,vacy)=0
 110	  continue
	  time(i+1)=atime
c-------------------------------------
c	Calculating sro and lro parameters
c-------------------------------------
	  do 140 ii=1,ns
	    do 150 jj=1,ns
	      if (site(ii,jj).lt.0) sro(i+1)=sro(i+1)+real(nn(ii,jj))
 150	    continue
 140	  continue
          sro(i+1)=sro(i+1)/(real(ns)**2)    ! P_B*P_BA
cNon-standard
c         sro(i+1)=sro(i+1)- RP_BA           ! diff. normalization, c=0 OK
cWCdef    
          sro(i+1)= 1 - sro(i+1)/RP_BA       ! 1 - P_BA/P_A 
	  do 160 jj=1,ns
	    do 170 ii=1,ns,2
	      if (mod(jj,2).eq.0) then
	        if (site(ii+1,jj).eq.1) lro(i+1)=lro(i+1) +1.
	      else
	        if (site(ii  ,jj).eq.1) lro(i+1)=lro(i+1) +1.
	      end if
 170	    continue
 160	  continue
	  lro(i+1)=abs(2.*lro(i+1)/real(natot)-1.)
 100	continue
ccccccccccccccccccccccc
c	Writing output
ccccccccccccccccccccccc
	write(6,*) '    temperature                bconc'
	write(6,*) tempk,bconc
	write(6,*) '       time	               sro                     lro'
	do 200 i=1,nout+1
	  write(6,*) time(i),sro(i),lro(i)
 200	continue
c----------------------
c	Storing the last configuration
c-----------------------
	open (unit=7,file='ordera.res', status='new')
	open (unit=8,file='orderv.res', status='new')
	open (unit=9,file='orderb.res', status='new')
	open (unit=10,file='order.res', status='new')
	do 300 i=1,ns
	  do 310 j=1,ns
	    if (site(i,j).eq.1) then
	      write(7,*) i,j
	    else if (site(i,j).eq.0) then
	      write(8,*) i,j
	    else
	      write(9,*) i,j
	    end if
 310	  continue
 300	continue
          do 350 j=1,ns
              write(10,'(1x,50i2)') (site(i,j), i=1,ns)
 350      continue
	close(7)
	close(8)
	close(9)
	stop
	end	
	
	double precision function ranno(ix)
	integer*4 a,p,ix,b15,b16,xhi,xalo,leftlo,fhi,k
	real*8 inverse
	parameter (inverse=4.65661287d-10)
	parameter (a=16807)
	parameter (b15=32768)
	parameter (b16=65536)
	parameter (p=2147483647)
	xhi=ix/b16
	xalo=(ix-xhi*b16)*a
	leftlo=xalo/b16
	fhi=xhi*a+leftlo
	k=fhi/b15
	ix=(((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
	  if (ix.lt.0) ix=ix+p
	ranno=dfloat(ix)*inverse
	continue
	return
	end  		




