/* * ising.c -- Very simple numerical Ising model simulation * * Written by Stefan Kesselheim 05/2009 for lecture tutorials. * Comments and marginal help by Florian Ruehle and Florian Dommert. * * ABOUT: * ------ * * Simulation of ISING model in a quadratic 2d area of varable length with * external magnetic field switched off (H=0). * The Ising model is an area of spins which point either up or down. * Nearest neighbour interaction is assumed (i.e. each spin has 4 neighbours), * periodic boundary conditions. * * LICENSE: * -------- * * This file may be distributed and/or modified under the terms of the * GNU General Public License version 2 as published by the Free Software * Foundation. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ #include /* Allocate memory */ #include /* Input and output */ #include /* GSL Random number generator */ #include /* For some calculations */ #include /* Read in time to seed RNG */ #include /* Read in process ID to seed RNG */ #define spin_type int /* Spin variable (basically 0 or 1) */ #define size_type int /* Size of the field */ #define magn_type long int /* Variable type for Magnetization */ /*************************** * Variables * ***************************/ gsl_rng* RNG; /* Random number generator */ size_type length_x; /* Field size in x direction */ size_type length_y; /* Field size in y direction */ int no_of_sweeps; /* How many time steps? */ int dump_every_nth_state; /* How often to dump a PPM */ int hot_cold_flag; /* Hot (1) or cold (0) start */ double beta; /* Variable for temperature adjustments */ double mag_sum; double mag_sq_sum; magn_type magnetization; /* Magnetization */ spin_type* field; /* Field of spins. Size adjusted at runtime. */ /*************************** * Functions * ***************************/ spin_type getSpin (size_type x, size_type y); /* Returns the current value of the spin at position (x,y) */ void setSpin (size_type x, size_type y, spin_type spin); /* Sets the spin at position (x,y) to value "spin" */ int flip (size_type x, size_type y); /* Flips the spin at position (x,y) to opposite and returns energy difference */ void create_field (unsigned short int mode); /* Allocates memory for the spin field, warm (mode=1) or cold (mode=0) start */ void delete_field (); /* Frees the memory used by the field */ void print_field (char* filename); /* Prints out the current state to a PPM file called "filename" */ int evaluate_mc_step (int change_of_energy); /* Decides whether to flip (returns 1) a spin or keep it (returns 0), depending on the change of energy */ void mc_sweep (); /* Sweeps through the whole field and flips spin if accepted. */ void run (); /* Initialize and run the simulation */ void initialize_RNG (); /* Seeds the RNG in the beginning */ int main (int argc, char** argv); /* Reads parameters from command line, assigns variables and calls run() */ /******************************** * Implementation * ********************************/ spin_type getSpin(size_type x, size_type y) { return field[ ((x + 10*length_x)%length_x)*length_x + ((y + 10*length_y)%length_y) ]; /* convert (x,y) to linear index first, then read out field-array */ } void setSpin(size_type x, size_type y, spin_type spin) { field[ ((x + 10*length_x)%length_x)*length_x + ((y + 10*length_y)%length_y) ] = spin; /* convert (x,y) to linear index first, then assign field-array */ } int flip (size_type x, size_type y) { int spin_sum_neighbourhood = 0; /* Variable for the sum of the four neighbour spins */ size_type no_of_parallel_neighbourhood_spins; /* Variable: How many are parallel to current spin? */ int change_of_energy; /* Variable: How much would the energy change in case of a flip? */ spin_sum_neighbourhood += getSpin(x+1,y); spin_sum_neighbourhood += getSpin(x-1,y); spin_sum_neighbourhood += getSpin(x,y+1); spin_sum_neighbourhood += getSpin(x,y-1); /* Add all spins in the neighborhood*/ /********* Do the flip: *********/ if (getSpin(x,y) == 0) { no_of_parallel_neighbourhood_spins = 4-spin_sum_neighbourhood; setSpin(x,y,1); magnetization += 2; } /* If spin is 0, flip to 1 and calculate No. of formerly parallel spins */ else { no_of_parallel_neighbourhood_spins = spin_sum_neighbourhood; setSpin(x,y,0); magnetization -= 2; } /* If spin is 1, flip to 0 and calculate No. of formerly parallel spins */ change_of_energy = 4*no_of_parallel_neighbourhood_spins -8; return change_of_energy; /* The amount of Delta(Energy) is returned */ } void create_field ( unsigned short int mode ) { spin_type spin; size_type i,j; field = malloc( length_x * length_y * sizeof(spin_type)); /* allocate the memory for the field at runtime */ magnetization = -length_x * length_y ; /* Arrange Magnetization around 0 */ if (mode == 0) for (i = 0; i < length_x; i++) for (j = 0; j < length_y; j++) { field[i*length_x+j] = 0; /* Cold start: All spins are 0 */ } if (mode == 1) for (i = 0; i < length_x; i++) for (j = 0; j < length_x; j++) { spin = (int) (2*gsl_rng_uniform (RNG)); field[i*length_x+j] = spin; /* Warm start: Spins are assigned randomly */ magnetization += 2*spin; } } void delete_field () { free(field); } void print_field(char* filename) { int i,j; FILE* file; file = fopen(filename, "w"); /* Open new file to write */ fprintf(file, "P1\n#\n%i %i\n", length_x, length_y); /* Write PPM header */ for ( i = 0; i < length_x; i++) { for ( j = 0; j < length_x; j++) if (getSpin(i,j) == 1) fprintf( file, "1 "); else fprintf( file, "0 "); fprintf( file, "\n"); /* Write all field data into file */ } fclose(file); } int evaluate_mc_step (int change_of_energy) { /* * Insert the Metropolis algorithm here!!! * A useful variable is "change_of_energy", which is sent to the function. * You will need a random number from "gsl_rng_uniform(RNG)". * * The function should return 1 if the flip is accepted * and 0 if it is declined. */ if (change_of_energy < 0) return 1; else /*missing part*/ return 1; else return 0; } void mc_sweep () { int i,j; for (i = 0; i < length_x; i++) for (j = 0; j < length_y; j++) { /* For all cells */ if (!evaluate_mc_step(flip(i,j))) flip(i,j); /* flip(i,j) returns the energy for the check, but also already flips the spin */ /* So, if evaluate_mc_step() returns 0, do another flip back to original state */ } } void run () { /* First, start random number generator and create field of wanted size: */ initialize_RNG(); create_field(hot_cold_flag); int i; char name[13]; /* for file name */ FILE* mag_file; mag_file = fopen("magnetization.dat","w"); printf("m, , \n"); for (i = 0; i < no_of_sweeps; i++) { /* for number of steps */ mc_sweep(); /* first do a monte carlo sweep */ fprintf(mag_file, "%f\n", (double) magnetization / length_x / length_y); if (dump_every_nth_state != 0 && i % dump_every_nth_state == 0) { /* If it's time to dump a picture */ sprintf(name, "dump%04i.ppm", i); /* assign file name */ print_field(name); /* and write file */ } mag_sum += (double) magnetization/ length_x / length_y; mag_sq_sum += (double) magnetization*magnetization/ length_x / length_y / length_x / length_y; } printf("%f %f %f\n", (double)magnetization/length_x/length_y,(double) mag_sum/no_of_sweeps, mag_sq_sum/no_of_sweeps); fclose(mag_file); /* Clean up: */ delete_field(); gsl_rng_free(RNG); } void initialize_RNG () { RNG = gsl_rng_alloc (gsl_rng_taus); /* allocate memory */ double seed = time (NULL) * getpid(); /* read in seed number from date/time and process ID */ gsl_rng_set (RNG, seed); /* seed the random number generator */ } int main ( int argc, char** argv ) { if (argc != 6) /* wrong number of parameters! */ exit(fprintf(stderr,"Usage: ising \n")); /* Print out how to use program */ else { /* Assign all variables from given parameters: */ beta = 1./atof(argv[1]); no_of_sweeps = atoi(argv[2]); length_x = atoi(argv[3]); length_y = length_x; hot_cold_flag = atoi(argv[4]); dump_every_nth_state = atoi(argv[5]); run(); /* Run the simulation */ } return 0; /* If the program doesn't crash, exit without error :o) */ }