// Monte Carlo simulation with the
// Metropolis scheme for spinodal decomposition

import java.lang.*;
import java.util.*;
import java.io.*;
import java.awt.*;  // uses the abstract windowing toolkit


public class Vacancy {
  // change the following globals to suit your fancy:
  static int maxSize = 200;      // width of window in pixels
  static int dim = 100;          // lattice size; can be any divisor of maxSize
  static int squareWidth = maxSize / dim;
  static Graphics g;  // the frame's graphics object, for drawing into
  static int frameTop, frameLeft;  // platform-dependent width of frame's top and left edges
  static Color aColor = new Color(00,160,60);
  static Color bColor = new Color(255,245,220);
  static Color vColor = new Color(150, 222, 250);

  static int seed;	
  static int ic; //counter for the current iteration
  static final int nsize = 10000;
  static final int nskip = 10000; 
  static final int pskip = 100; //replot lattice after pskip iterations
  static final int ntotal = nsize*nskip; //total number of iterations
  static final int neq = 0;	//# of iterations for equilibriation
  static final int lx = dim; //dimensions of lattice
  static final int ly = dim; //dimensions of lattice
  static int iaccept = 0;  // counter for calculating the acceptance rate of attempted jumps
  static int at; 
  static int a[][] = new int[lx][ly]; // array that stores the nature of species in the lattice
  static Random r = new Random();
  static final int rmax = dim/2+1; //maximum value of correlation length
  static final double jaa = 1;  
  static final double jab = -1; 
  static final double Teq = 1.0; //for equilibriation at temperature Teq before final simulations
  static double T1 = 1.5;		
  static double T;
  static double r1;  
  static double fV = 0.001;		//fraction of vacant sites
  static double fA = 0.5; 		// fraction of A atoms
  static double fB = 1.0 - fV - fA; //fraction of B atoms
  
  public static void main(String argv[]) {
  // initialize window and graphics:
    Frame gWin = new Frame("Vacancy Model");
    gWin.setLocation(50,50);  // screen coordinates of top left corner
    gWin.setResizable(false);
    gWin.setVisible(true);  // show it!
    Insets theInsets = gWin.getInsets();
    gWin.setSize(maxSize+theInsets.left+theInsets.right,maxSize+theInsets.top+theInsets.bottom);
    frameTop = theInsets.top;
    frameLeft = theInsets.left;
    long resumeTime = System.currentTimeMillis() + 1000;
    do {} while (System.currentTimeMillis() < resumeTime);  // wait for frame to get initialized
    g = gWin.getGraphics();  // stick graphics object into a global so we can draw from anywhere!

  	
  	FileOutputStream out;
    PrintStream p;

  	// Initiate the seed from the current time
    GregorianCalendar t = new GregorianCalendar();
    int t1 = t.get(Calendar.SECOND);
    int t2 = t.get(Calendar.MINUTE);
    int t3 = t.get(Calendar.HOUR_OF_DAY);
    int t4 = t.get(Calendar.DAY_OF_MONTH);
    int t5 = t.get(Calendar.MONTH) + 1;
    int t6 = t.get(Calendar.YEAR);
	seed = t6+70*(t5+12*(t4+31*(t3+23*(t2+59*t1))));

    if((seed%2) == 0) seed = seed-1;
    
  	initial(a);			//initializes the lattice 
  	
  	System.out.println("\n" + "Initial lattice: ");	
	
    ///Open file for writing	
	
	try{
       out = new FileOutputStream("corVacancy.dat",false);
       p = new PrintStream(out);             
       p.print("\n");
       p.close();
    }
	catch(Exception e){
         System.err.println("Error in printing to file");
    }
    
	///////////////////////////////////////////////////////

  	double accept = 0;
  	
  	for (int i=0; i<neq; ++i){   //for equilibriating at Temperature Teq before quenching to lower temperature
  		iaccept = 0;
  		T = Teq;
   		exchange();  //exchange atom and vacancy positions between neighbours
  		accept = i*accept + iaccept*100.0/(lx*ly);
  		accept /= i+1;
  	}
  
    System.out.println("\n" + "After " + neq + " exchanges: ");
    
	System.out.println("Accept rate = " + accept + "%");
	
    double cr[] = new double[rmax+1];    
    
    for(int ic=0; ic<=ntotal; ++ic) {
    	iaccept = 0;
    	T = T1;  
        exchange(); //exchange atom and vacancy positions between neighbours
           		
   		accept = (neq+ic)*accept + iaccept*100.0/(lx*ly);
  		accept /= neq+ic+1;
  		
        if(ic%nskip == 0){
      		System.out.println("\n" + "ic = " + ic);  
			cr = correl();
			//append correlation function values for this iteration
			try{
		       out = new FileOutputStream("corVacancy.dat",true);
		       p = new PrintStream(out);             
		       p.print("\n");
		       p.println("ic = " + ic);
		       for(int r=0; r<=rmax; ++r){
		       	   p.println(r + " " + cr[r]);
		       } 
		       p.println();   
		       p.close();
		     }
		     catch(Exception e){
		           System.err.println("Error in printing to file");
		     }

      		 System.out.println("Accept rate = " + accept + "%");
        }        
		//replot the lattice 
		if(ic%pskip==0){
	    	for(int i=0; i<lx; ++i)
	        	for(int j=0; j<ly; ++j)
        			colorSquare(i,j);     		
    	} 
    }

  }

/////////////////////////////////////////////////////////  
//Initialize the lattice
  public static void initial(int a[][]){ 
  	
  	for(int i=0; i<lx; ++i)
  		for(int j=0; j<ly; ++j){
  			r1 = ranf();
  			if(r1<=fA)  a[i][j] = 1; // A atoms
  			else if(r1>=(fA+fV))  a[i][j] = -1;  // B atoms
  			else a[i][j] = 0; //Vacant site  	
  		}
  	
  	at = 0;			
  	int nab = 0;	//# of occupied sites
  	for(int i=0; i<lx; ++i)
  		for(int j=0; j<ly; ++j){
  			at += a[i][j];
  			nab += Math.abs(a[i][j]);
  			colorSquare(i,j);
  		}  
  		
  	int nV = lx*ly - nab;	//# of vacant sites
  	int cnt = (int) (fV*lx*ly) - nV;
  	
  	System.out.println(at + " more atoms of A than B");
    System.out.println("# of vacant sites = " + nV + " out of " + lx*ly + " lattice sites" ); //check
    //for small lattice size, there should be at least 1 vacancy.
 }
	
///////////////////////////////////////////////////////////	
	

//Diffusion
 public static void exchange(){
 int swp = 0;
 int i2=0, j2=0;
 
 double delH = 0;
 double kB = 1;
 double p; 
 int xaa=0, xab=0, yaa=0, yab=0;
 int daa=0, dab=0;
 
 iaccept = 0;
 //Exchange of vacancy with atoms 
 for(int i=0; i<lx; ++i)
 	for(int j=0; j<ly; ++j){
 		i = (int) (lx*ranf());
 		j = (int) (ly*ranf());
 		if(a[i][j]==0){		//if vacancy on site
 		//	select the neighbour (i2,j2) for exchange
 			int t1 = (int) (4.*ranf());
 			if(t1==0){
 				i2 = (lx+i-1)%lx;
 				j2 = j;
 			}
 			if(t1==1){
 				i2 = i;
 				j2 = (ly+j-1)%ly;
 			} 			 
 			if(t1==2){
 				i2 = (i+1)%lx;
 				j2 = j;
 			} 			 
 			if(t1==3){
 				i2 = i;
 				j2 = (j+1)%ly;
 			}
 			if(a[i2][j2]!=0){
 			// # of AA or BB bonds before exchange
	 			xaa =  Math.max(a[i2][j2]*a[(lx+i2-1)%lx][j2],0) + Math.max(a[i2][j2]*a[(i2+1)%lx][j2],0); 	    
		        xaa +=  Math.max(a[i2][j2]*a[i2][(ly+j2-1)%ly],0) + Math.max(a[i2][j2]*a[i2][(j2+1)%ly],0);
		    // # of AB bonds before exchange    
		        xab =  -Math.min(a[i2][j2]*a[(lx+i2-1)%lx][j2],0) - Math.min(a[i2][j2]*a[(i2+1)%lx][j2],0); 	    
		        xab +=  -Math.min(a[i2][j2]*a[i2][(ly+j2-1)%ly],0) - Math.min(a[i2][j2]*a[i2][(j2+1)%ly],0);
		    // # of AA or BB bonds after exchange    
		        yaa =  Math.max(a[i2][j2]*a[(lx+i-1)%lx][j],0) + Math.max(a[i2][j2]*a[(i+1)%lx][j],0); 	    
		        yaa +=  Math.max(a[i2][j2]*a[i][(ly+j-1)%ly],0) + Math.max(a[i2][j2]*a[i][(j+1)%ly],0) - 1;
		    // # of AA or BB bonds after exchange    
		        yab =  -Math.min(a[i2][j2]*a[(lx+i-1)%lx][j],0) - Math.min(a[i2][j2]*a[(i+1)%lx][j],0); 	    
		        yab +=  -Math.min(a[i2][j2]*a[i][(ly+j-1)%ly],0) - Math.min(a[i2][j2]*a[i][(j+1)%ly],0); 
	        
	            daa = yaa - xaa;  //change in A atoms
		        dab = yab - xab;  //change in B atoms 
		        delH = -(daa*jaa + dab*jab); //change in energy

		        p = Math.exp(-delH/(kB*T));
		        if(p>=ranf()){					//accept if delH<0 else,
		        	a[i][j] = a[i2][j2];		// except with a probability
		        	a[i2][j2] = 0;  
		        	++iaccept;   
		        } 				
 			} 		
 		} 	    
    }
 }
//////////////////////////////////////////////////////////////////////////////// 
  
 
 // Method to generate a uniform random number in [0,1]
// following x(i+1)=a*x(i) mod c with a=pow(7,5) and
// c=pow(2,31)-1.  Here the seed is a global variable.

  public static double ranf() {
    final int a = 16807, c = 2147483647, q = 127773,
      r = 2836;
    final double cd = c;
    int h = seed/q;
    int l = seed%q;
    int t = a*l-r*h;
    if (t > 0) seed = t;
    else seed = c+t;
    return seed/cd;
  }
 
////////////////////////////////////////////////////////////////////////////////
//method to calculate the correlation function

  public static double[] correl() {
 	double cr[] = new double[rmax+1]; 
 	double temp = 0.0;
	   
 	for(int r=0; r<=rmax; ++r){ 	
 	    cr[r] = 0.0;
	    for(int i=0; i<lx; ++i)
	    	for(int j=0; j<ly; ++j){
	    		if(r==0) cr[0] += a[i][j]*a[i][j];
	    		if(r!=0){
	    			temp = a[i][(j+r)%ly] + a[(i+r)%lx][j];
	    			cr[r] += a[i][j]*temp/2.0;	    			
	    		}   		 
	    	}	    
	}
	
	for(int r=0; r<=rmax; ++r){
		cr[r] /= lx*ly;
	}		
	
    return cr;
 }
 
////////////////////////////////////////////////////////////////////////////////
 // given a lattice site, color that square according to the current s value:
    private static void colorSquare(int i, int j) {
        if (a[i][j] == 1) g.setColor(aColor);
        if (a[i][j] == -1) g.setColor(bColor);
        if (a[i][j] == 0) g.setColor(vColor);
        g.fillRect(squareWidth*i + frameLeft, squareWidth*j + frameTop, squareWidth, squareWidth);
    }   
}