#define PROGNAME "mc_Ising_2D" #define VERSION "0.4" #define DATE "15.01.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): Wolf cluster updates */ #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 Wolf 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 rec(int mag[NMAX][NMAX], int Inn[4], int Jnn[4], int L, int *acc, int i, int j, double prob, int mag0) /* recursion subroutine for constructing the Wolf cluster */ { int k, Inew, Jnew,is,js; (*acc)++; for(k=0;k<4;k++){ // loop over nearest neighbors 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; /* for(is = 0; is < L; is++){ */ /* for(js = 0; js < L; js++) */ /* printf("%d ",cluster[is][js]); */ /* printf("\n"); */ /* } */ /* printf ("%d %d %d %d %d %d %d\n", i,j, Inew, Jnew, mag[i][j], mag[Inew][Jnew], cluster[Inew][Jnew]); */ if (cluster[Inew][Jnew]==0){ /* printf ("mag0*mag[Inew][Jnew]: %d\n", mag0*mag[Inew][Jnew]); */ if ((mag0*mag[Inew][Jnew]==1) && (drand48() < prob )) { cluster[Inew][Jnew]=1; mag[Inew][Jnew]=-mag[Inew][Jnew]; rec(mag,Inn, Jnn, L, acc,Inew,Jnew, prob, mag0); } } } return; } void mc_Wolf(int mag[NMAX][NMAX], int Inn[4], int Jnn[4], int L, double prob, double *time) { int i, j, acc; int mag0; /* initialize cluster */ for (i=0;i<L;i++) for (j=0;j<L;j++){ cluster[i][j]=0; /* printf("cluster: cluster[%d][%d]=%d\n",i,j,cluster[i][j]); */ } /* select seed site for cluster */ i = (int) (drand48()*L); j = (int) (drand48()*L); mag0=mag[i][j]; mag[i][j]=-mag[i][j]; cluster[i][j]=1; // outputcluster(L); /* printf("%d %d %d %lf\n",i,j,mag0, prob); */ acc=0; rec(mag, Inn, Jnn, L, &acc,i,j, prob, mag0); *time+=1.0*acc/(L*L); 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)olf cluster\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 Wolf */ 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 Wolf 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: %2.1f\n", 100.0*accept/move); } else { // Wolf updates for(i=0;i<warm;i++) mc_Wolf(mag, Inn, Jnn, L, prob, &simtime); for(i=0;i<sweeps;i++){ mc_Wolf(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 }