/*******************************************************************************

   Cluster.c

This program is to simulate 3D Ising model using Wolff's cluster algorith[1],
for which oLsqy a single cluster is constructed and updated at each sweep.
A single cluster is constructed and updated at each sweep. From the starting
configuration, we choose a site at random and construct a cluster around it
by bonding together neighboring sites with the appropriate probabilities.
All sites in this cluster are then given the same new spinvalue, producing
the new configuration, which is obviously far less correlated with the initial
configuration than the result of a single Metropolis[2] update.
Wolff's method is probably the best sequential cluster algorithm.

 [1] U. Wolff, Phys. Rev. Lett. 62,361(1989).
 [2] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller,
 j. Chem. Phys. 21, 1087(1953).

  Written 10/28/00 by Liping Huang

*******************************************************************************
*******************************************************************************

   List of Integers
 
   L                                   The size of the L-cube lattice.
   Lsq                                 Number of spinper layer.
   NSpin                               Total Number of spinin 3D.   
   NBin                                Number of bins into which each 
                                       run is broken into.
   NPass                               Number of passes 		
   NPassBin                            Number of passes per bin
   NPassEq                             Number of equilibriation passes
   NTemp                               Number of Temperatures.


   List of Interger Array

   SetI[NSpin]                         Set index indicating which cluster the site is in                  
   SetS[NSpin]                         Set size
   Spin[NSpin]                         The Spin at i-th Lattice site
   NbrLst[NSpin][6]                    Neighbor table

   List of Floating Variables

   dE                                  Energy difference
   E                                   Energy
   Esq                                 Energy squared
   M                                   Magnetic Moment
   Msq                                 Moment squared
   T                                   Temperature
 
   List of Floating Arrays

   EBin[NBin]                          The Energy
   MBin[NBin]                          The Magnitization
   Cv[NBins]                           The Heat Capacity  
   X[NBins]                            The Susceptibility Chi
  
   List of Fucntion

   void bond(int, int);
   int search(int);
********************************************************************************/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"RDMG"               /* Random number generator from numerical recipes.*/

#define L 10              /* Some of these will be read from a file */
#define Lsq L*L
#define NSpin L*L*L
#define NBin 10
#define Nkstep 500
#define NPass 5000

int SetI[NSpin];                     
int SetS[NSpin];   
       
main()
{
 
  int Spin[NSpin], NbrLst[NSpin][6]; 
  
  int i,j,k,m,i0, nbr;
  int conf, intM,intE;
  int ispin,ibin,ipass,iK;

  double E,Esq,M,Msq,Cv,Cvsq,X,Xsq,E_bin,Esq_bin,M_bin,Msq_bin,dE;
  double DeltaE,DeltaM,DeltaCv,DeltaX,NPassBin,NPassEq,K,dK,fNBin;
  double CvBin[NBin],XBin[NBin];
  double EBin[NBin],MBin[NBin];
  double prob, Pr[3][13];
 
  void bond(int, int);
  int search(int);
  
  FILE *ifp;


    ifp=fopen("C-M-bin10.dat","a");

    /* Output parameters */  
	
    fprintf(ifp, " L is: %d, NBin is: %d, NPass is: %d,  Nkstep is: %d\n",L,NBin,NPass,Nkstep);

    /* Write column headings */

    fprintf(ifp, " K      Energy    DeltaE    Magnet    DeltaM    SpHt      DeltaCv   Suscep    DeltaX\n");
    fclose(ifp);
  

    /* Define some useful parameters */

    NPassBin=NPass/NBin;
    NPassEq=NPass/10;
  

    /* Initialize a spin-up system */
   
   for(i = 0; i < NSpin; i++)
       Spin[i]= 1;
         
 
/* Construct the Neighbour Table */

   for( i = 0; i < L; i++)
   for( j = 0; j < L; j++)
   for( k = 0; k < L; k++)
    {
      m = i *Lsq + j*L + k;
      NbrLst[m][0] = i*Lsq + j*L +(k +1)%L;
      NbrLst[m][1] = i*Lsq + j*L +(k + L -1)%L;
      NbrLst[m][2] = i*Lsq + ((j+1)%L)*L + k;
      NbrLst[m][3] = i*Lsq + ((j + L - 1 )%L)*L + k;
      NbrLst[m][4] = ((i + 1)%L)*Lsq + j*L + k;
      NbrLst[m][5] = ((i + L-1)%L)*Lsq + j*L + k;
     }



  K=0.15; dK=0.0003;

  
/* Carry out the simulation at a fixed temperature starting with T */
   K-=dK;
   

  for(iK=1;iK<=Nkstep;iK++)
  {
    K+=dK;
    /* Calculate the probability array for Metropolis algorithm */

     for(j=0;j<=2;j+=2)
      {
        for(k=0;k<=12;k++) /* There are 13 different possibilities */
        { 
	      dE=(float)(2*(j-1)*(k-6));  
          Pr[j][k]=exp(-K*dE);     /* Hence Pr[-1]=Pr[0]; Pr[+1]=Pr[2] */
        }
      }



    /* Do some equilibration passes (NPassEq=NPass/10) */
    
   if(K<0.25)  //Eualibrate with Wolf Cluster algorithm.
   {
       printf("Wolf Cluster running at the temperature: %f \n", K);
       /*Calculate the probability to build bonds between neighbors.*/ 
       prob = exp(-2.0*K);
   

/* Do some equilibration passes (NPassEq=NPass/10) */

    for(ipass=1;ipass<=NPassEq;ipass++)
     {
     	for(i=0;i<NSpin;i++)
            {
              SetI[i]=i;
              SetS[i]=1;
             }
           
         
         for(i=0;i<NSpin;i++)
            {  m = (int)(i/Lsq);
               j = (int)((i-m*Lsq)/L);
               k = (int)(i-m*Lsq-j*L);
          if(k!=(L-1))
		       {
                if(Spin[i]== Spin[NbrLst[i][0]] && rand2()>prob)
                   bond(i,NbrLst[i][0]);
                }
	      if(j!=(L-1))
               { 
                if(Spin[i]== Spin[NbrLst[i][2]] && rand2()>prob)
                   bond(i,NbrLst[i][2]);
               }
          if(m!=(L-1))
              {
                if(Spin[i]== Spin[NbrLst[i][4]] && rand2()>prob)
                   bond(i,NbrLst[i][4]);
              }
           }

          i0 = (int)(rand2()*NSpin);
          i0=search(i0);


          for(i=0;i<NSpin;i++)
            {
              if(SetI[i]==i0)
              Spin[i] = -Spin[i];
            }
          }
        }
  else //Eualibrate with Metropolis algorithm.
   {
     printf("Metropolis running at the temperature: %f \n", K);

    
    for(ipass=1;ipass<=NPassEq;ipass++)
      {
        for(i=0;i<NSpin;i++)
         {
          conf=6;           /* Such that -6 maps to 0 of the array */
	      for(nbr=0;nbr<=5;nbr++) /* ie. conf = 0 is Spin = -6 */
	      conf+=Spin[NbrLst[i][nbr]]; /* This ranges from -6 to 6 */
	      if(rand2()<Pr[Spin[i]+1][conf])
	      Spin[i]=-Spin[i];
		}
       }
    }

/* Break the run into NBin bins Loop through the bins */
   for(ibin=1;ibin<=NBin;ibin++)
    {
      /* Initialization */
      E_bin=0.0;         /* Total Energy in the bin */
      Esq_bin=0.0;       /* Energy squared */       
      M_bin=0.0;         /* Magnetic Momemnt in the bin */
      Msq_bin=0.0;       /* Moment squared */
      CvBin[ibin]=0.0;   /* Specific heat */
      XBin[ibin]=0.0;
      EBin[ibin]=0.0;
      MBin[ibin]=0.0;


/* Now the passes which count are done for thSetI particular j-th bin */
     
     for(ipass=1;ipass<=NPassBin;ipass++)
	 {   /*update the Spins. */
        for(i=0;i<NSpin;i++)
         {
          conf=6;           /* Such that -6 maps to 0 of the array */
	      for(nbr=0;nbr<=5;nbr++) /* ie. conf = 0 is Spin = -6 */
	      conf+=Spin[NbrLst[i][nbr]]; /* This ranges from -6 to 6 */
	      if(rand2()<Pr[Spin[i]+1][conf])
	      Spin[i]=-Spin[i];
	 }
       
/* Calculate M and E after each pass */

	   intM=0;
	   intE=0;

	 for(i=0; i<NSpin; i++)
	  {
	    intM +=Spin[i];   /* This is the M for one pass */
        intE -=Spin[i]*(Spin[NbrLst[i][0]]+Spin[NbrLst[i][2]]+Spin[NbrLst[i][4]]);
	  }


/* Accumulate the Energy, Magnetic moment and their squares */
	
	  E_bin+=(float)(intE);
	  M_bin+=(float)(intM);
	  Esq_bin+=(float)(pow(intE,2));
	  Msq_bin+=(float)(pow(intM,2));
    }  //end of pass loop.



      /* Compute for the bin just finished, the average of the above calculated quantities */
      E_bin=E_bin/(float)(NPassBin);
      M_bin=M_bin/(float)(NPassBin);
      Esq_bin=Esq_bin/(float)(NPassBin);
      Msq_bin=Msq_bin/(float)(NPassBin);
      /* Determine for this bin, the heat capacity per Spin C, susceptibility X, energy per spin */
      
	  CvBin[ibin]=(Esq_bin-pow(E_bin,2))*K*K/(float)(NSpin);
      XBin[ibin]=(Msq_bin-pow(M_bin,2))*K/(float)(NSpin);
      EBin[ibin]=E_bin/(float)(NSpin);
      MBin[ibin]=M_bin/(float)(NSpin);
   } //end of bin loop.


    /* All passes at the current temperature are complete; compute averages
       over the entire run and determine uncertainties 

       Initialisation of the sums */

    E=0.0;   Cv=0.0;    M=0.0;    X=0.0;
    Esq=0.0; Cvsq=0.0;  Msq=0.0;  Xsq=0.0;

    /* Accumulate the results from different bins */
  for(ibin=1;ibin<=NBin;ibin++)
   {
      Cv+=CvBin[ibin];
      X+=XBin[ibin];
      E+=EBin[ibin];
      M+=MBin[ibin];
      Cvsq+=pow(CvBin[ibin],2);
      Xsq+=pow(XBin[ibin],2);
      Esq+=pow(EBin[ibin],2);
      Msq+=pow(MBin[ibin],2);
   }

    /* Normalise the sums over the bins */
    fNBin=(float)(NBin);
   
    Cv/=fNBin;
    Cvsq/=fNBin;
    X/=fNBin;
    Xsq/=fNBin;
    E/=fNBin;
    Esq/=fNBin;
    M/=fNBin;
    Msq/=fNBin;

    /* Calculate the uncertainties */
    
    DeltaCv=sqrt((Cvsq-Cv*Cv)/(float)(NBin-1));
    DeltaX=sqrt((Xsq-X*X)/(float)(NBin-1));
    DeltaE=sqrt((Esq-E*E)/(float)(NBin-1));
    DeltaM=sqrt((Msq-M*M)/(float)(NBin-1));
    
    /* Output the results */
    
    ifp=fopen("C-M-bin10.dat","a");
    fprintf(ifp,"%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n",K,E,DeltaE,M,DeltaM,Cv,DeltaCv,X,DeltaX);
    fclose(ifp);

 } //end of Temperature loop.


}  //end of main().

void bond(int i, int j)
{
    	
  while(SetI[i]!=i)
  i=SetI[i];
  while(SetI[j]!=j)
  j=SetI[j];
  if(i==j) return;
  if(SetS[i]<SetS[j])
     {
       SetS[j]=SetS[j]+SetS[i];
       SetI[i]=j;
      }
  else
   {
     SetS[i]=SetS[i]+SetS[j];
     SetI[j]=i;
   }
}

int search(int i)
{
  while(SetI[i]!=i)
  i=SetI[i];
  return i;
}