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

                                                                             
   PROGRAMME: 3D Ising Model ( Metropolis Algorithm )
									


   Written                                                      11/11/00
   Modified                                                  
     1.Analysis by Histogram Technique K=K0                     12/03/00
     2.exp(-DeltaK*E) used instead                              12/06/00

   Dyutiman Das
   Dept. of Physics
   University of Illinois, Urbana,IL

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

   List of Integers

   conf                                Configuration in the probability
                                       matrix.
   L                                   The size of the L-cube lattice.
   Lsq                                 Number of spins per layer.
   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
   NSpin                               Number of spins (+/-).           
   NTemp                               Number of Temperatures.


   List of Integer Arrays

   NbrSpin[i][j]                       The spin of the j-th neighbour
                                       of the i-th spin.
   Spin[NSpin]                         The Spin at i-th Lattice site



   List of Floating Variables

   dE                                  Energy difference
   E                                   Energy
   Esq                                 Energy squared
   M                                   Magnetic Moment
   Msq                                 Moment squared
   T                                   Temperature
 




   List of Floating Arrays

   Cv[NBins]                           The Heat Capacity  
   Pr[SpinValue][Conf]                 The probability Array
   X[NBins]                            The Susceptibility Chi

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

									
#include<stdio.h>
#include<math.h>
#include"RDMG"

#define L 20              /* Some of these will be read from a file */
#define Lsq L*L
#define NSpin L*L*L
#define NE 6*NSpin
#define NBin 10
#define Nkstep 500
#define NPass 500000
#define dK 0.0001
#define Ksim 0.224


main(){

  int i,j,k,NLayer,conf,intE,intM,iE,iM,kstep;
  int ipass,ispin,nbr,NPassEq,ibin;
  int Spin[NSpin],NbrLst[NSpin][6];
  int H[NE+1],m[NE+1],msq[NE+1],absm[NE+1],m4[NE+1];

  double dE,fE,fM,DeltaK,norm,Mnorm,K,HexpE;
  double fNPass,fNSpin,Pr[3][13],M,E,Msq,Esq;
  double Mabs,M4,dlnmsqdk,dUdk,dUdk_Max;
  double absmE,msqE,absm_Max,msq_Max,dlnabsmdk;

  double fNBin,Cv_Max,Xt_Max,Mabs_Max,Msq_Max;
  double spheat,suscep,dma,dmb;
  double avKc_Cv,avKcsq_Cv,DKc_Cv;
  double avKc_Xt,avKcsq_Xt,DKc_Xt;
  double avKc_Mabs,avKcsq_Mabs,DKc_Mabs;
  double avKc_Msq,avKcsq_Msq,DKc_Msq;
  double av_m2,av_m,D_m,av_msq2,av_msq,D_msq;


  double Kc_Cv[NBin],Kc_Xt[NBin],Kc_Mabs[NBin],Kc_Msq[NBin];
  double Max_m[NBin],Max_msq[NBin];


  FILE *ifp;



  /* Get the input parameters from the input file */
  /* IN FACT THESE ARE IN THE ABOVE #define COMMAND */

  fNBin=(float)(NBin);
  fNPass=(float)(NPass);
  fNSpin=(float)(NSpin);
  Mnorm=fNSpin*fNPass;
  NPassEq=NPass/10;

  /* Initialise the lattice to be Spin UP */
  for(ispin=0;ispin<NSpin;ispin++)
    Spin[ispin]=1;

  /* Construct the Neighbour Table */
  for(ispin=0;ispin<NSpin;ispin++){
    NbrLst[ispin][0]=ispin+1;   NbrLst[ispin][1]=ispin-1;   /* East & West */
    NbrLst[ispin][2]=ispin+L;   NbrLst[ispin][3]=ispin-L;   /* North & South */
    NbrLst[ispin][4]=ispin+Lsq; NbrLst[ispin][5]=ispin-Lsq; /* Up & Down */    
  }
  /* Fix the boundary spins */
  for(ispin=0;ispin<Lsq;ispin++){
    NbrLst[ispin][5]=(L-1)*Lsq+ispin;  /* (L-1) layers, each with Lsq lattice
                                          points, plus the current point from 
					  the beginning of that layer; this 
					  fixes the periodic BC of the bottom 
					  layer, 5 = Down */
    NbrLst[(L-1)*Lsq+ispin][4]=ispin;  /* This fixes the BC of the top layer */
  }
  for(ispin=0;ispin<NSpin;ispin+=L){
    NbrLst[ispin][1]=ispin+L-1;          /* Left face, 1 = West, */
    NbrLst[ispin+L-1][0]=ispin;          /* Right face */
  }
  for(j=0;j<=NSpin-Lsq+1;j+=Lsq){
    for(k=0;k<L;k++){
      ispin=j+k;
      NbrLst[ispin][3]=ispin+L*(L-1);    /* South face */
      NbrLst[ispin+L*(L-1)][2]=ispin;    /* North face */
    }
  }






  /* Carry out the simulation at a fixed temperature Ksim */

  /* Calculate the probability array */
  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(-Ksim*dE);     /* Hence Pr[-1]=Pr[0]; Pr[+1]=Pr[2] */
    }
  }
  /* Do some equilibration passes (NPassEq=NPass/10) */
  for(ipass=1;ipass<=NPassEq;ipass++){
    /* Run sequentially through the lattice, flipping
	 spins as appropriate */
    for(ispin=0;ispin<NSpin;ispin++){
      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[ispin][nbr]]; /* This ranges from -6 to 6 */
      if(rand2()<Pr[Spin[ispin]+1][conf])
	Spin[ispin]=-Spin[ispin];
    }
  }
  

  /* DISTRIBUTE THE RUN INTO NBIN BINS */

  for(ibin=0;ibin<NBin;ibin++){
    /* First initialise the histogram  and the m averages */
    for(i=0;i<=NE;i++){
      H[i]=0;m[i]=0;msq[i]=0;m4[i]=0;
    }
    
    /* Initailise the max values and the K arrays */
    Cv_Max=0.0;      Xt_Max=0.0;      Mabs_Max=0.0;      Msq_Max=0.0;
    Kc_Cv[ibin]=0.0; Kc_Xt[ibin]=0.0; Kc_Mabs[ibin]=0.0; Kc_Msq[ibin]=0.0;

    for(ipass=1;ipass<=NPass/NBin;ipass++){
      /* Run sequentially through the lattice flipping spins
	 as appropriate */
      for(ispin=0;ispin<NSpin;ispin++){
	conf=6;           
	for(nbr=0;nbr<=5;nbr++)
	  conf+=Spin[NbrLst[ispin][nbr]]; /* This ranges from -6 to 6 */
	if(rand2()<Pr[Spin[ispin]+1][conf])
	  Spin[ispin]=-Spin[ispin];
      }
      /* Calculate M and E after each pass */
      intM=0.0;
      intE=0.0;
      for(ispin=0;ispin<NSpin;ispin++){
	intM+=(Spin[ispin]);   /* This is the M for one pass */
	for(nbr=0;nbr<=4;nbr+=2)
	  intE-=(Spin[ispin]*Spin[NbrLst[ispin][nbr]]); 
      }
      iE=intE+3*NSpin;
      H[iE]+=1;
      m[iE]+=intM;
      absm[iE]+=abs(intM);
      msq[iE]+=intM*intM;
      m4[iE]+=pow(intM,4);
    }
    
    /* The histogram has been constructed,

       the analysis of the data (acquired at one temperature) is 
       carried out by the histogram technique. */

    for(kstep=0;kstep<=Nkstep;kstep++){ /* kstep, dK,*/
      DeltaK=(float)(kstep-Nkstep/2)*dK; 
      K=Ksim+DeltaK; 
      /* Calculate the norm for the probability */
      norm= 0.0;
      for(i=0;i<=NE;i++){
	iE=i-3*NSpin;
	norm+=(float)(H[i])*exp(-DeltaK*(float)(iE));
      }

      /* Calculate the averages of the basic quantities */
      E=0.0;Esq=0.0;M=0.0;Msq=0.0;Mabs=0.0;M4=0.0;
      for(i=0;i<=NE;i++){
	iE=i-3*NSpin; 
	HexpE=(float)(H[i])*exp(-DeltaK*(float)(iE));
	E+=(float)(iE)*HexpE;
	Esq+=(float)(iE*iE)*HexpE;
	M+=(float)(m[i])*HexpE/Mnorm;
	Msq+=(float)(msq[i])*HexpE/Mnorm;
	Mabs+=(float)(absm[i])*HexpE/Mnorm;
	M4+=(float)(m4[i])*HexpE/Mnorm;
      }
  
      E/=norm;Esq/=norm;    /* Normalise the average */
      M/=norm;Msq/=norm;    /*    This is the histogram normalisation */
      Mabs/=norm;M4/=norm;  /*    not over passes. */

      /* HERE WE HAVE TO CALCULATE ALL THE REQUIRED THERMODYNAMIC
	 QUANTITIES */
      absmE=0.0;msqE=0.0;
      /* Calculate m^2E_av, |m|E_av */
      for(i=0;i<=NE;i++){
	iE=i-3*NSpin;
	HexpE=(float)(H[i])*exp(-DeltaK*(float)(iE));
	absmE+=(float)(absm[i]*iE)*HexpE/Mnorm; 
	msqE+=(float)(msq[i]*iE)*HexpE/Mnorm;
      }
      absmE/=norm; msqE/=norm; /* Normalise */


      /* Specific Heat and Susceptibility */
      spheat=(Esq-pow(E,2))*K*K/fNSpin;
      suscep=(Msq-pow(M,2))*K/fNSpin;
      dma=(absmE/Mabs)-E;  /* dln|m|/dK */
      dmb=(msqE/Msq)-E;    /* dlnm^2/dK */

      if(spheat>Cv_Max){
	Cv_Max=spheat; Kc_Cv[ibin]=K;
      }
      if(suscep>Xt_Max){
	Xt_Max=suscep; Kc_Xt[ibin]=K;
      }
      if(Mabs>Mabs_Max){
	Mabs_Max=Mabs;    Kc_Mabs[ibin]=K;
      }
      if(Msq>Msq_Max){
	Msq_Max=Msq;   Kc_Msq[ibin]=K;
      }
      
      if(dma>absm_Max)
	absm_Max=dma;
      if(dmb>msq_Max)
	msq_Max=dmb;

    }                                    /* K loop ends */    
    Max_m[ibin]=absm_Max;
    Max_msq[ibin]=msq_Max;
  }                                      /* Bin loop ends */

  av_m=0.0;       av_msq=0.0;
  avKc_Cv=0.0;    avKcsq_Cv=0.0;
  avKc_Xt=0.0;    avKcsq_Xt=0.0;
  avKc_Mabs=0.0;  avKcsq_Mabs=0.0;
  avKc_Msq=0.0;   avKcsq_Msq=0.0;

  for(ibin=0;ibin<NBin;ibin++){
    av_m+=Max_m[ibin];      av_m2+=pow(Max_m[ibin],2);
    av_msq+=Max_msq[ibin];  av_msq2+=pow(Max_msq[ibin],2);

    avKc_Cv+=Kc_Cv[ibin];avKcsq_Cv+=pow(Kc_Cv[ibin],2);
    avKc_Xt+=Kc_Xt[ibin];avKcsq_Xt+=pow(Kc_Xt[ibin],2);
    avKc_Mabs+=Kc_Mabs[ibin];avKcsq_Mabs+=pow(Kc_Mabs[ibin],2);
    avKc_Msq+=Kc_Msq[ibin];avKcsq_Msq+=pow(Kc_Msq[ibin],2);
  }

  av_m/=fNBin;    av_m2/=fNBin;
  av_msq/=fNBin;  av_msq2/=fNBin;

  avKc_Cv/=fNBin; avKcsq_Cv/=fNBin;
  avKc_Xt/=fNBin; avKcsq_Xt/=fNBin;
  avKc_Mabs/=fNBin; avKcsq_Mabs/=fNBin;
  avKc_Msq/=fNBin; avKcsq_Msq/=fNBin;


  D_m=sqrt(fNBin*abs(av_m2-pow(av_m,2))/(fNBin-1.0));
  D_msq=sqrt(fNBin*abs(av_msq2-pow(av_msq,2))/(fNBin-1.0));

  DKc_Cv=sqrt(fNBin*(avKcsq_Cv-pow(avKc_Cv,2))/(fNBin-1.0));
  DKc_Xt=sqrt(fNBin*(avKcsq_Xt-pow(avKc_Xt,2))/(fNBin-1.0));
  DKc_Mabs=sqrt(fNBin*(avKcsq_Mabs-pow(avKc_Mabs,2))/(fNBin-1.0));
  DKc_Msq=sqrt(fNBin*(avKcsq_Msq-pow(avKc_Msq,2))/(fNBin-1.0));

  ifp=fopen("Kc.dat","a");
  fprintf(ifp,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",L,avKc_Cv,DKc_Cv,avKc_Xt,DKc_Xt,avKc_Mabs,DKc_Mabs,avKc_Msq,DKc_Msq);
  fclose(ifp);

  ifp=fopen("nu.dat","a");
  fprintf(ifp,"%d\t%f\t%f\t%f\t%f\n",L,av_m,D_m,av_msq,D_msq);   
  fclose(ifp); 

}