/* sa.c This is search for minimum energy spin glass configurations, using simulated annealing. On a DEC Alpha, be sure to compile this program with the -fptm sui and -g flags, otherwise the floating point exception from the addition of a series of random numbers of zero mean and unit variance will cause it to halt quickly ! J. Watlington, 11/26/95 */ #include #include #include #define REPORT_PERIOD 0x7F /* num of time steps between prints */ #define SPIN_GLASS_SIZE 100 /* number of spin glass lattice points */ #define EQUAL_PROB 0.5 #define MAX_GLASS (SPIN_GLASS_SIZE - 1) /* The following variables represent the system being modeled. */ double glass_spin[ SPIN_GLASS_SIZE ]; /* The spin states */ double glass_J[ SPIN_GLASS_SIZE ]; /* The external stresses */ long rseed = -1; /* random generator state */ /* Quitting conditions */ double epsilon = 1.0; #define MAX_STALE_PASSES 200 #define MAX_ITERATIONS 400000 /* The following Numerical Recipes in C routines should be linked in. */ extern float gasdev( long *idnum ); extern float ran1( long *idnum ); /* The following routines are declared locally */ int accept_transition( double delta_E, double beta ); double calc_total_energy( void ); void display_glass( int time, int show_J ); void main( int argc, char **argv ) { double alpha; double beta = 0.0; double total_E; double delta_E; double minimum_E = 0.0; int time = 0; int i, x; int quitting_time = 0; /* The sole argument is a cooling rate, alpha */ if( argc != 2 ) { printf( "usage: sa \n" ); exit( -1 ); } sscanf( argv[1], "%lf", &alpha ); /* Initialize the spin glass */ for( i = 0; i < SPIN_GLASS_SIZE; i++ ) { glass_J[ i ] = (double)gasdev( &rseed ); minimum_E = minimum_E - fabs( glass_J[ i ] ); } for( i = 0; i < SPIN_GLASS_SIZE; i++ ) { if( ran1( &rseed ) > EQUAL_PROB ) glass_spin[ i ] = 1.0; else glass_spin[ i ] = -1.0; } total_E = calc_total_energy(); /* Tell the user about the starting conditions */ printf( "Starting energy: %lf\n", total_E ); printf( "Looking for minimum energy %lf (to within %lf)\n", minimum_E, epsilon ); minimum_E = minimum_E + epsilon; display_glass( time, 0 ); /* Now start iterating */ do { /* Pick a point in the spin glass, and calculate the energy change contributed by a spin flip */ x = (int)(ran1( &rseed ) * (float)MAX_GLASS ); if( x == MAX_GLASS ) delta_E = -((glass_J[x-1] * glass_spin[x-1] * -glass_spin[x]) + (glass_J[x] * -glass_spin[x] * glass_spin[0])) +((glass_J[x-1] * glass_spin[x-1] * glass_spin[x]) + (glass_J[x] * glass_spin[x] * glass_spin[0])); else if( x == 0 ) delta_E = -((glass_J[MAX_GLASS] * glass_spin[MAX_GLASS] * -glass_spin[x]) + (glass_J[x] * -glass_spin[x] * glass_spin[x+1])) +((glass_J[MAX_GLASS] * glass_spin[MAX_GLASS] * glass_spin[x]) + (glass_J[x] * glass_spin[x] * glass_spin[x+1])); else delta_E = -((glass_J[x-1] * glass_spin[x-1] * -glass_spin[x]) + (glass_J[x] * -glass_spin[x] * glass_spin[x+1])) +((glass_J[x-1] * glass_spin[x-1] * glass_spin[x]) + (glass_J[x] * glass_spin[x] * glass_spin[x+1])); if( accept_transition( delta_E, beta ) ) { glass_spin[ x ] = -glass_spin[ x ]; total_E += delta_E; quitting_time = 0; } else quitting_time++; beta += alpha; /* increment beta (beta = alpha * time) */ if( !(time & REPORT_PERIOD) ) { printf( "t: %4d energy: %lf\n", time, total_E ); #ifdef FAILSAFE_CHECK delta_E = calc_total_energy(); if( fabs( total_E - delta_E ) > 0.001 ) { printf( "internal error: energy loss (%f)!\n", total_E-delta_E ); total_E = delta_E; } #endif } } while( (total_E > minimum_E) && (time++ < MAX_ITERATIONS) && (quitting_time < MAX_STALE_PASSES) ); if( quitting_time >= MAX_STALE_PASSES ) printf( "Glass frozen in non-optimal state !\n" ); printf( "Ending energy: %lf at time %d\n", total_E, time ); display_glass( time, 1 ); } int accept_transition( double delta_E, double beta ) { float threshold; if( delta_E < 0.0 ) { return( 1 ); } threshold = (float)exp( -beta * delta_E ); if( ran1( &rseed ) < threshold ) return( 1 ); else return( 0 ); } double calc_total_energy( void ) { int x; double total; for( x = 0; x < SPIN_GLASS_SIZE-1; x++ ) { total += glass_J[x] * glass_spin[x] * glass_spin[ x+1 ]; } /* Handle periodic boundary condition */ x = SPIN_GLASS_SIZE - 1; total += glass_J[x] * glass_spin[x] * glass_spin[ 1 ]; return( -total ); } #define NUM_COLS 25 void display_glass( int time, int show_J ) { int index = 0; int row; int col; static int num_rows = SPIN_GLASS_SIZE / NUM_COLS; printf( "At time %d:", time ); if( show_J ) { printf( "\n**** J" ); for( row = 0; row < num_rows; row++ ) { printf( "\n%2d>", row ); for( col = 0; col < NUM_COLS; col++ ) { printf( " %2.1f", glass_J[ index++ ] ); } } printf( "\n**** S" ); index = 0; } for( row = 0; row < num_rows; row++ ) { printf( "\n%2d>", row ); for( col = 0; col < NUM_COLS; col++ ) { printf( " %2.0f", glass_spin[ index++ ] ); } } printf( "\n" ); }