#define PROGNAME "mc_Ising_2D"
#define VERSION "0.4t"
#define DATE  "10.12.2007"
#define AUTHOR "Nils Bluemer"

/* Metropolis Monte Carlo for 2D-Ising Model with pbc */
/* Code adapted from http://gold.cchem.berkeley.edu/~acpan/220A/2DIsing.c */
 /* NEW (0.3): code clean-up */
 /* NEW (0.4): Wolff cluster updates */
/* template based on 0.4b */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*************************
 * Constant Declarations *
 ************************/

#define NMAX 50            /* Lattice size (i.e. NMAX x NMAX) */

#define J -1.0             /* Ising lattice coupling constant
			    J < 0: Ferromagnetic
			    J > 0: Antiferromagnetic */
#define H 0.0              /* Ising lattice field strength */

 /* global variables */
  int cluster[NMAX][NMAX];       /* flag for Wolff cluster updates */

void error(char error_text[])
/* standard error handler */
{
    fprintf(stderr,"\n%s run-time error\n",PROGNAME);
    fprintf(stderr,"--%s--\n",error_text);
    fprintf(stderr,"for general help use option -h\n");
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}


void random_init(int mag[NMAX][NMAX], int L)
{
  int i,j;

  for(i = 0; i < L; i++)
    for(j = 0; j < L; j++){
      if(drand48() > 0.5)
	mag[i][j] = 1;
      else
	mag[i][j] = -1;
    }
}

void up_init(int mag[NMAX][NMAX], int L)
{
  int i,j;

  for(i = 0; i < L; i++)
    for(j = 0; j < L; j++){
	mag[i][j] = 1;
    }
}

void outputmag(int mag[NMAX][NMAX], int L)
{
  int i,j;

  printf("\n");

  for(i = 0; i < L; i++){
    for(j = 0; j < L; j++){
      if(mag[i][j] == -1)
	printf("* ");
      else
	printf("  ");
    }
    printf("\n");
  }
}

void outputcluster(int L)
{
  int i,j;

  printf("\n");

  for(i = 0; i < L; i++){
    for(j = 0; j < L; j++){
      if(cluster[i][j] == -1)
	printf("* ");
      else
	printf("  ");
    }
    printf("\n");
  }
}

double mag_av(int mag[NMAX][NMAX], int L)
/* computes average magnetization */
{
  int i,j;
  double x;

  x=0.0;

  for(i = 0; i < L; i++)
    for(j = 0; j < L; j++)
      x+=1.0*mag[i][j];
  x=x/(L*L);
  return(x);
}

double compute_E(int mag[NMAX][NMAX], int Inn[4], int Jnn[4], int L)
/* computes total energy for given spin configuration */
{
  int i, j, Inew, Jnew, k;
  double Energy;
  Energy = 0.0;

  for(i=0;i<L;i++)
    for(j=0;j<L;j++)
      {
	/* Loop over nearest neighbors */
	for(k=0;k<=1;k++)
	  {
	    Inew = i + Inn[k];
	    Jnew = j + Jnn[k];

	    /* Check periodic boundary conditions */
	    if(Inew < 0)
	      Inew = L-1;
	    else if(Inew >= L)
	      Inew = 0;

	    if(Jnew < 0)
	      Jnew = L-1;
	    else if(Jnew >= L)
	      Jnew = 0;

	    /* Update the energy */
	    Energy += J * mag[i][j] * mag[Inew][Jnew];
	  }
	/*Calculate the contribution from the field H */
	Energy += H*mag[i][j];
      }
  return(Energy);
}

void mc_single_spin (int mag[NMAX][NMAX], int Inn[4], int Jnn[4], int L,
		     int *accept, int *move, double *Energy, double beta)
/* performs one sweep of single spin-flips (random choice of spins) */
{
  double Etemp, deltaE;
  int n, k, i, j, Inew, Jnew;

  for(n=0;n<L*L;n++)
    {
      i = (int) (drand48()*L);
      j = (int) (drand48()*L);

      Etemp = 0.0;

      for(k=0;k<4;k++){
	Inew = i + Inn[k];
	Jnew = j + Jnn[k];

	if(Inew < 0)
	  Inew = L-1;
	else if(Inew >= L)
	  Inew = 0;

	if(Jnew < 0)
	  Jnew = L-1;
	else if(Jnew >= L)
	  Jnew = 0;

	Etemp += -J*mag[i][j]*mag[Inew][Jnew];
      }

      deltaE = 2.0*Etemp - 2*H*mag[i][j];

      if(deltaE < 0){
	  mag[i][j] = -mag[i][j];
	  *Energy += deltaE;
	  (*accept)++;
	}
      else if(drand48()<exp(-beta*deltaE)){
	  mag[i][j] = -mag[i][j];
	  *Energy += deltaE;
	  (*accept)++;
	}
      (*move)++;
    }
}

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

void mc_Wolff(int mag[NMAX][NMAX], int Inn[4],
	     int Jnn[4], int L, double prob, double *time)
{
  return;
}

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

void printhelp ()
{
 printf("**********************************************************\n");
 printf("%s: Metropolis Monte Carlo for 2D Ising model\n",PROGNAME);
 printf("Version: %s, %s by %s\n",VERSION,DATE,AUTHOR);
 printf("based on http://gold.cchem.berkeley.edu/~acpan/220A/2DIsing.c\n");
 printf("options: -T temperature (really: k_B T/|J|)\n");
 printf("         -L linear dimension of lattice (L<=%d)\n",NMAX);
 printf("         -w# number of initialization sweeps\n");
 printf("         -n# number of measurement sweeps\n");
 printf("         -c print configurations\n");
 printf("         -u$ update: (s)ingle-spin or (W)olff cluster\n");
 printf("         -i initialization with all spins up (not random)\n");
 printf("         -h this help\n");
}

/***********************************************************/
int main(int argc, char *argv[])
{

  /*************************
   * Variable Declarations *
   ************************/

  int mag[NMAX][NMAX];           /* 2D Ising Lattice */
  int i;                         /* Loop counters */
  double Energy;                 /* Total lattice energy */
  int Inn[4] = {1,0,-1,0};       /* Nearest neighbor array I */
  int Jnn[4] = {0,1,0,-1};       /* Nearest neighbor array J */
  int accept = 0;                /* Number of accepted moves */
  int move = 0;                  /* Number of moves total */

  char c;
  double T;                      /* temperature (in units of J/k_B) */
  double beta;
  int sweeps;                    /* number of measurement sweeps */
  int warm;                      /* number of warm-up sweeps */
  int print_conf;                /* flag for printing configurations */
  int L;                         /* lattice dimension */

  char update, init;                   /* flag for update: single-spin or Wolff */
  double prob, simtime=0.0;


  /***************************
   * Initialization          *
   ***************************/

  T=1.5;
  sweeps=500;
  print_conf=0;
  L=NMAX;
  update='s';
  init='r';

  while (--argc > 0 && (*++argv)[0] == '-')
    while (c= *++argv[0])
	   switch (c) {
	   case 'n':
 	     sscanf(++argv[0],"%d\n",&sweeps);
             break;
	   case 'w':
 	     sscanf(++argv[0],"%d\n",&warm);
             break;
	   case 'L':
 	     sscanf(++argv[0],"%d\n",&L);
             if (L>NMAX) error ("L too large");
             break;
           case 'T':
 	     sscanf(++argv[0],"%lf\n",&T);
             break;
           case 'u':
 	     sscanf(++argv[0],"%c\n",&update);
             break;
	   case 'i':
	     init='u';
	     break;
           case 'c':
             print_conf=1;
             break;
	   case 'h':
	     printhelp();
	     exit(0);
             break;
/*  	   default:  */
/*  	     error("No valid choice");  */
	   }

  /* compute derived variables */
  beta=1.0/T;
  prob=(1.0-exp(2.0*beta*J)); // for Wolff cluster updates

  /* print header */
  printf("# L=%d, T=%f, w=%d, n=%d, print_conf=%d\n", L, T, warm, sweeps, print_conf);

  /* Seed the random number generator */
  srand48((unsigned int) time(NULL) - (time(NULL)/100000)*100000);

  /* Choose initial spin configuration */
  if (init=='r')
    random_init(mag,L);
  else
    up_init(mag, L);

  /* Determine the energy initially */
  Energy=compute_E(mag,Inn,Jnn,L);

  if (update=='s'){ // single-spin flips
    for(i=0;i<warm;i++)
      mc_single_spin (mag, Inn, Jnn, L, &accept, &move, &Energy, beta);
    for(i=0;i<sweeps;i++){
      mc_single_spin (mag, Inn, Jnn, L, &accept, &move, &Energy, beta);
      printf ("%9.6lf  %9.6lf\n", mag_av(mag,L),Energy/(L*L));
    }

    /*Ratio of accepted spin flips compared to total attempted flips*/
    printf("# Acceptance rate: %5.3f\n", 1.0*accept/move);
  } else { // Wolff updates
    for(i=0;i<warm;i++)
      mc_Wolff(mag, Inn, Jnn, L, prob, &simtime);
    for(i=0;i<sweeps;i++){
      mc_Wolff(mag, Inn, Jnn, L, prob, &simtime);
      Energy=compute_E(mag,Inn,Jnn,L);
      printf ("%9.6lf  %9.6lf  %9.6lf\n", mag_av(mag,L),Energy/(L*L), simtime);
    }
  }

  if (print_conf>0) outputmag(mag, L); // print final configuration

}