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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Structure of DATA file:  order.data
c
c10000000   20    1   ! #vac jumps, #sub-blocks, init conf.(1=ord,3=rdm)
c  380 	  0.50   1.0  ! temp.in K, concentration, Ballistic Rate (0-10)
c  0.0    -0.05  0.0  ! Vaa, Vab, Vbb for first n.n. shell
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	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,jx,jy,jn
	integer*4 iseed,nstep
	integer*4 g1,sitet,incr,count
	real*8 gball,g0,x1,y1,RP_BA
	real*8 vaa,vab,vbb,temp,tempk,bconc
	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 count /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	g0    : Scaled rate  (range =0-10.0, w/ 0.01 a low athermal rate)
c	gball : Rate for Ballistic Jumps (here forced swap of atoms)
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,g0
	read (5,*) vaa,vab,vbb
	temp=tempk/11604.9
	close (5)
ccccccccccccccccccccccc
c	Initialization
ccccccccccccccccccccccc
c       Athermal Ballistic jump frequencies (scale set so that g0 input 
c                      is O(1) and gball increases as system increases)
        gball= g0*real(nssq)/9000000.0  
c--------------------------------
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 directly 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)
          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.0 - sro(1)/RP_BA             ! 1 - P_BA/P_A
c
	  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        Athermal Ballistic jump frequencies = gball
c        i.e., not governed by thermal process.
c---------------------------------------
c--------------------------------------
c	Choosing the jump to be performed
c--------------------------------------
c  write(6,*) ' w-s',w(1),w(2),w(3),w(4)
          dtime=1/(w(1)+w(2)+w(3)+w(4)+gball)
	  wcumul=0.
              g1=0
	  ran=ranno(iseed)/dtime
	  do 120 k=1,4
	    wcumul=wcumul+w(k)
            if (ran.le.wcumul) goto 125
 120	  continue
            g1=1
 125	  continue
	  atime=atime+dtime
c--------------------------------------
c	Updating neighborhoods and positions
c--------------------------------------
       if (g1.gt.0 .and. gball.gt.0.0) then     ! athermal jump
c
        count=count+1
          jx=int(ranno(iseed)*real(ns)) + 1  ! arb. site
          jy=int(ranno(iseed)*real(ns)) + 1
          jn=int(4.0*ranno(iseed) + 0.999999999)      ! arb. neighbor
           x1=modl(jx+nvec(1,jn))
           y1=modl(jy+nvec(2,jn))
c
c Do something if two sites are different
          if(site(jx,jy) .ne. site(x1,y1)) then
c
c Is first or second site an A, if so add one in proper place? 
            incr=0
            if (site(jx,jy).eq.1) then
              incr= 1
            elseif(site(x1,y1).eq.1) then
              incr=-1
            endif
 
               do 1301 kk=1,4
                 xx=modl(jx+nvec(1,kk))
                 yy=modl(jy+nvec(2,kk))
                 nn(xx,yy)=nn(xx,yy) - incr
 1301          continue
               do 1351 kk=1,4
                 xx=modl(x1+nvec(1,kk))
                 yy=modl(y1+nvec(2,kk))
                 nn(xx,yy)=nn(xx,yy) + incr
 1351          continue

              if(site(x1,y1).eq.0) then
                  vacx=jx
                  vacy=jy
              end if
              if(site(jx,jy).eq.0) then
                  vacx=x1
                  vacy=y1
              end if
            
c If different, exchange atoms in site array
          sitet=site(jx,jy)
          site(jx,jy)=site(x1,y1)
          site(x1,y1)=sitet
    
          endif
c
       else                  ! original thermal jump
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
c 
       endif
c
 110	  continue
c-------------------------------------
c       Finished NSTEPS  update NSTEP time
c-------------------------------------
	  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.0 - sro(i+1)/RP_BA         ! 1 - P_BA/P_A 
c
	  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,*) ' Athermal Jumps ', count,' out of ',nstep*nout,' Jumps'
	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)
	close(10)
	close(11)
	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