/* * xising * * Xwindow program to simulate the 2-d Ising model. The underlying * approach uses many "demons" circulating around the lattice trying * to flip spins. * * The latest version of this program can be found at * http://thy.phy.bnl.gov/www/xtoys/xtoys.html * * Michael Creutz creutz@bnl.gov * November 1999, version 2.11: corrects problems on alphas pointed out by Kari Rummukainen * * compiling: cc -O -o xising -L/usr/X11R6/lib xising.c -lm -lX11 * * On suns: * * cc -O -I/usr/openwin/include -L/usr/openwin/lib newxising.c -lm -lX11 * * If you find some machine supporting X that this does not work on, * please let me know. * */ /* the usual includes */ #include #include #include #include #include #include #include #include #include /* * initial lattice dimensions: ncols will be truncated to a multiple * of WORDSIZE nrows will be changed if a multiple of 29 since I let * the demons hop 29 sites to get a bit of randomness */ int nrows = 256, ncols = 256; long **field = NULL; /* to store the system; it will be * dynamically allocated to a [ncols] * by [nrows/WORDSIZE] array */ /* some convenient definitions */ #define WORDSIZE (CHAR_BIT*sizeof(long)) #define LEFTBIT (~LONG_MAX) #define nextrow (row+1-nrows*(row>=nrows-1)) /* nshift is here to allow for different endianness */ #define nextcol (col+nshift*(1-hwords*(col==topcol))) #define previousrow (row-1+nrows*(row<1)) #define previouscol (col-nshift*(1-hwords*(col==bottomcol))) char stringbuffer[100]; static char *progname; char bits[256]; /* for counting bits */ /* hwords will store the number of long words in a single row */ int hwords, volume; /* bitcount counts the set bits in a word; call setup() before using */ #define bitcount8(i) (bits[(i)&255]) #define bitcount16(i) (bitcount8(i)+bitcount8((i)>>8)) /* check if long int is 32 bits */ #if ((0x7FFFFFFF)==(LONG_MAX)) /* 32 bits */ #define bitcount(i) (bitcount16(i)+bitcount16((i)>>16)) /* note, current linux glibc has a bug in mrand48, so it only gives 16+1 random bits. When that is fixed, this can go back to mrand48() */ #define RND (mrand48()^lrand48()) #else /* assume 64 bits */ #define bitcount32(i) (bitcount16(i)+bitcount16((i)>>16)) #define bitcount(i) (bitcount32(i)+bitcount32((i)>>32)) #define RND (mrand48()^(mrand48()<<32)) #endif /* things to monitor temperature, etc. */ double beta; /* the inverse temperature */ long energycount = 0, volumecount = 0; int heater = 0, cooler = 0; int canonical = 1; int boundary = 0, algorithm = 0, paused = 0, iteration = 0; /* variables to allow for peculiar bit orderings on some machines */ int nshift, topcol, bottomcol; /* for local algorithm; use WORDSIZE two-bit demons */ long demon0, demon1, checkermask; long *work0 = NULL, *work1 = NULL; /* * bitprob represents the probability for the corresponding demon bit * to be set multiplied by 1<<16 and rounded to an integer; used to * speed up demon refreshing */ long bitprob0, bitprob1; /* stuff for the cluster algorithm; use three demon bits */ #define DEMONBITS 3 int ndemons; long **cluster = NULL; /* to store cluster shape */ long **sadx = NULL, **sady = NULL; /* labels sad demons */ int **unchecked = NULL; /* to reduce redundant calculations */ int boxtop, boxbottom; /* will enclose cluster */ long **demon = NULL; /* for demons in cluster algorithm */ long **newdemon = NULL; int *activerow = NULL; long demonindex = 0; /* for shuffling demons */ long clustergrowth = 0; int direction = 1; /* for vertical sweeps through * cluster */ /* things for regulating the updating speed */ int speed = 0; /* updating delay proportional to * speed */ int delay = 0; /* counter for implementing speed * control */ struct timeval timeout; /* timer for speed control */ /* various window stuff */ Display *display; int screen; Window window, quitbutton, pausebutton, algobutton, heatwindow, boundwindow, speedbutton, canonbutton, makebutton(); GC gc, gcpen, gcclear; XImage *spinimage = NULL, *clusterimage = NULL; XFontStruct *font = NULL; int font_height, font_width; unsigned int windowwidth, windowheight; XSizeHints size_hints; int depth, darkcolor, lightcolor, black, white, spinup, spindown; long event_mask; /* all my functions */ void refreshdemons(int step); /* refreshes the cluster demons, step * steps over some of them for only a * partial refresh */ void drawbutton(), openwindow(), setup(), makebuttons(), growcluster(), startcluster(), localupdate(), repaint(), setalg(), setheater(), setboundary(), fixthermometer(), check(), leftshift(), rightshift(), allocarray(), heatbath(); /* text for the buttons */ char *boundtext[4] = {"periodic", "antiperiodic", "up boundary", "down"}; char *heatertext[3] = {"run free", "heat", "cool"}; char *algotext[2] = {"local", "cluster"}; /* general layout dimensions */ #define LEFTOFFSET 64 #define TOPOFFSET 124 #define HEATERHEIGHT 54 #define HEATERWIDTH 86 #define BOUNDHEIGHT 72 #define BOUNDWIDTH 116 #define BUTTONWIDTH 48 #define BUTTONHEIGHT 18 #define TTOP 128 #define THGT (nrows*2) #define TBOTTOM (TTOP+THGT) #define TLEFT 24 #define TWDTH 5 int main(argc, argv) int argc; char **argv; { unsigned int width, height, i; XEvent report; progname = argv[0]; openwindow(argc, argv); setup(); makebuttons(); /* loop forever, looking for events */ while (1) { if (clustergrowth) /* finish cluster regardless of * algorithm */ while (clustergrowth) { growcluster(); } else if ((0 == XPending(display)) && (0 == paused)) { /* no events waiting */ if (delay) { /* don't update yet */ delay--; /* * this use of select() seems a kludge to me; why can't * usleep() be more standard? */ timeout.tv_sec = 0; timeout.tv_usec = 50000;/* .05 sec per delay unit */ select(0, NULL, NULL, NULL, &timeout); } else { delay = speed; /* reset delay counter */ if (1 == algorithm) { startcluster(); } else localupdate(); } } if (paused | XPending(display)) { XNextEvent(display, &report); /* find out what X wants done */ switch (report.type) { case Expose: /* the window has been exposed and * needs redrawing */ if ((report.xexpose.window) != window) break; /* cuts down flashing, but you might * remove this line if things aren't * being redrawn */ if (report.xexpose.count != 0) break; /* more in queue, wait for them */ repaint(); break; case ConfigureNotify: /* the window has been resized */ width = report.xconfigure.width; height = report.xconfigure.height; if ((width < size_hints.min_width) || (height < size_hints.min_height)) { fprintf(stderr, "%s: window too small to proceed.\n", progname); XUnloadFont(display, font->fid); /* I'm still too neat */ XFreeGC(display, gc); /* from my Amiga days */ XFreeGC(display, gcpen); XFreeGC(display, gcclear); XCloseDisplay(display); exit(1); } if ((width != windowwidth) || (height != windowheight)) { windowwidth = width; windowheight = height; /* new lattice dimensions */ ncols = width - 75; nrows = (height - 164 - 24) / 2; makebuttons(); } break; case ButtonPress: if (report.xbutton.window == quitbutton) { XUnloadFont(display, font->fid); XFreeGC(display, gc); XFreeGC(display, gcpen); XFreeGC(display, gcclear); XCloseDisplay(display); exit(1); } else if (report.xbutton.window == pausebutton) { paused = 1 - paused; drawbutton(pausebutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT, "pause", 1 - 2 * paused); } else if (report.xbutton.window == algobutton) { /* toggle algorithm */ algorithm = 1 - algorithm; /* 1 for cluster, 0 for local updating */ setalg(); for (i = 0; i < 2; i++) drawbutton(algobutton, 0, i * 18, 68, 18, algotext[i], 1 - 2 * (algorithm == i)); } else if (report.xbutton.window == heatwindow) { /* adjust heater */ setheater((report.xbutton.y * 3) / HEATERHEIGHT); } else if (report.xbutton.window == boundwindow) /* set boundary conditions */ { setboundary((report.xbutton.y * 4) / BOUNDHEIGHT); } else if (report.xbutton.window == canonbutton) { canonical = 1 - canonical; if (canonical) drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT, "canonical", -1); else drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT, "microcanonical", 1); } else if (report.xbutton.window == speedbutton) { /* * reset speed; speed=0 is the fastest, 10 the slowest */ speed = 10 - (11 * (report.xbutton.x - 1)) / (BOUNDWIDTH - 2); if (speed < 0) speed = 0; if (speed > 10) speed = 10; delay = speed; drawbutton(speedbutton, 0, 0, BOUNDWIDTH, 18, "speed", -1); drawbutton(speedbutton, 1 + ((10 - speed) * (BOUNDWIDTH - 2)) / 11, 1, (BOUNDWIDTH - 2) / 11, 16, "", 2); } else { /* do an update on a random mouse * press */ if (canonical && (report.xbutton.y > TOPOFFSET) && (report.xbutton.x < LEFTOFFSET)) { beta = 0.5 / (0.6 + 1.0 * (TBOTTOM - report.xbutton.y) / THGT); bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); fixthermometer(); } else { if (1 == algorithm) startcluster(); else localupdate(); } } break; default: break; } } } return 1; } void startcluster() { /* flip the old cluster and start a * new one */ long xedge, yedge; int row, col, bit; long currentword; long index0, index1; /* flip previous cluster and adjust demons on boundary */ for (row = 0; row < nrows; row++) { for (col = 0; col < hwords; col++) { index0 = (demonindex + 2 * (col + row * hwords)); index1 = (index0 + 1); while (index0 >= ndemons) index0 -= ndemons; /* faster than using % */ while (index1 >= ndemons) index1 -= ndemons; if (0 == unchecked[row][col]) { /* don't bother if unchecked */ unchecked[row][col] = 1;/* reset unchecked */ currentword = cluster[row][col]; /* find cluster edges */ xedge = currentword ^ ((currentword << 1) | (1 & (cluster[row][nextcol] >> (WORDSIZE - 1)))); yedge = currentword ^ cluster[nextrow][col]; /* update demons and field */ for (bit = 0; bit < DEMONBITS; bit++) { demon[bit][index0] = (newdemon[bit][index0] & xedge) | ((~xedge) & demon[bit][index0]); demon[bit][index1] = (newdemon[bit][index1] & yedge) | ((~yedge) & demon[bit][index1]); } field[row][col] ^= currentword; } /* end of if checked */ /* do any cooling or heating required */ /* * This doesn't work very well in the cluster algorithm; the * demons hold the heat for a while before dumping it, thus the * thermomenter doesn't read reliably. For major heating and * cooling it seems better to go to the local algorithm. */ if ((!canonical) && (drand48() < 0.05)) { if (heater) { demon[0][index0] ^= (RND & RND); demon[0][index1] ^= (RND & RND); } else if (cooler) { demon[0][index0] &= RND; demon[0][index1] &= RND; } else { /* rotate the demon's first bit to * help scramble */ demon[0][index0] = (demon[0][index0] << 1) | (1 & (demon[0][index0] >> (WORDSIZE - 1))); demon[0][index1] = (demon[0][index1] << 1) | (1 & (demon[0][index1] >> (WORDSIZE - 1))); } } /* accumulate counts */ currentword = demon[0][index0]; energycount += bitcount(currentword); currentword = demon[0][index1]; energycount += bitcount(currentword); volumecount += WORDSIZE; } /* end of col loop */ } /* end of row loop */ if (canonical) { if (heater || cooler) { beta += (cooler - heater) * .0005; bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); fixthermometer(); } } /* clear previous cluster */ for (row = 0; row < nrows; row++) { activerow[row] = 0; for (col = 0; col < hwords; col++) cluster[row][col] = 0; } /* randomize starting demon location */ demonindex = lrand48() % ndemons; /* pick random site to seed cluster */ row = lrand48() % nrows; col = lrand48() % hwords; clustergrowth = 1; cluster[row][col] = 1 << (lrand48() % WORDSIZE); activerow[row] = activerow[previousrow] = activerow[nextrow] = 1; /* set up bounding box for cluster */ boxbottom = boxtop = row; /* monitor temperature */ if (canonical) { refreshdemons(5); /* refresh 1/10th of the demons */ } else if (volumecount >= volume) { /* update beta and * thermometer */ beta = .5 * log((double) (-1. + volumecount / (0.5 * energycount))); bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); fixthermometer(); energycount = 0; volumecount = 0; } /* end of temperature monitoring */ /* copy lattice to window */ XPutImage(display, window, gc, spinimage, 0, 0, LEFTOFFSET, TOPOFFSET, ncols, nrows); /* display the new cluster */ XPutImage(display, window, gc, clusterimage, 0, 0, LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows); return; } /* end of startcluster */ void growcluster() /* this new version (July `97) is supposed to be more readable */ { int rowloop, row, col, nrow, ncol, newstuff; long edge, currentword, newsites; clustergrowth = 0; direction = -direction; /* change sweeping direction */ /* start where last activity was */ if (direction > 0) row = boxtop; else row = boxbottom; for (rowloop = 0; rowloop < 4 * nrows; rowloop++) { /* multiple passes seem to help */ if (activerow[row]) { activerow[row] = 0; nrow = nextrow; for (col = 0; col < hwords; col++) { ncol = nextcol; currentword = cluster[row][col]; /* look for new stuff if part of cluster is around */ newstuff = (0 != (currentword | cluster[nrow][col] | cluster[row][ncol])); while (newstuff) { newstuff = 0; /* make sure necessary demons calculated */ if (unchecked[row][col]) check(row, col); /* now put things from/to next column */ if (1 & (sadx[row][col])) { if (1 & (currentword ^ ((cluster[row][ncol]) >> (WORDSIZE - 1)))) { newstuff = 1; cluster[row][ncol] |= LEFTBIT; currentword |= 1; } } /* put things from/to next row */ edge = currentword ^ (cluster[nrow][col]); if ((newsites = edge & sady[row][col])) { /* single "=" intended */ currentword |= newsites; cluster[nrow][col] |= newsites; newstuff = 1; activerow[nrow] = 1; } /* find cluster edge within current word */ edge = (currentword ^ (currentword << 1)) & (~1); while ((newsites = edge & (sadx[row][col]))) { /* need to grow */ currentword |= newsites | (newsites >> 1); edge = (currentword ^ (currentword << 1)) & (~1); newstuff = 1; } if (newstuff) { cluster[row][col] = currentword; activerow[row] = activerow[previousrow] = 1; clustergrowth = 1; } } /* end of while newstuff */ } /* end of column loop */ } /* end of if activerow */ if (direction > 0) { /* alternate directions */ if (activerow[row]) /* back up */ boxtop = row = nextrow; else row = previousrow; } else { if (activerow[row]) /* back up */ boxbottom = row = previousrow; else row = nextrow; } if (speed) /* speed is not an issue; draw * progress */ XPutImage(display, window, gc, clusterimage, 0, row, LEFTOFFSET, TOPOFFSET + 28 + nrows + row, ncols, 1); } /* end of row loop */ /* display new cluster */ XPutImage(display, window, gc, clusterimage, 0, 0, LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows); return; } /* end of growcluster */ void check(row, col) int row, col; { /* calculate sadx and sady when * unchecked */ /* demons are sad if they cannot flip a bond */ long right, bottom, carry0, carry1; int bit; long index0, index1; right = field[row][col] ^ ((field[row][col] << 1) | (1 & (field[row][nextcol] >> (WORDSIZE - 1)))); bottom = field[row][col] ^ field[nextrow][col]; if (boundary) { /* correct for antiperiodic * boundaries */ if (row == nrows - 1) bottom ^= (-1); if (col == topcol) right ^= 1; } index0 = (demonindex + 2 * (col + row * hwords)); /* start of xdemons */ index1 = (index0 + 1); /* start of ydemons */ while (index0 >= ndemons) index0 -= ndemons; while (index1 >= ndemons) index1 -= ndemons; /* copy demons to newdemon */ for (bit = 0; bit < DEMONBITS; bit++) { newdemon[bit][index0] = demon[bit][index0]; newdemon[bit][index1] = demon[bit][index1]; } /* subtract 1 from newdemon */ carry0 = carry1 = (-1); for (bit = 0; bit < DEMONBITS; bit++) { newdemon[bit][index0] ^= carry0; newdemon[bit][index1] ^= carry1; carry0 &= newdemon[bit][index0]; carry1 &= newdemon[bit][index1]; } /* add two if neighbor antiparallel */ for (bit = 1; bit < DEMONBITS; bit++) { newdemon[bit][index0] ^= right; newdemon[bit][index1] ^= bottom; right &= (~newdemon[bit][index0]); bottom &= (~newdemon[bit][index1]); } carry0 ^= right; carry1 ^= bottom; /* make demons sad if overflow */ sadx[row][col] = carry0; sady[row][col] = carry1; unchecked[row][col] = 0; return; } void setup() { int i, j, count; /* initialize "bits" used for bitcounts */ for (i = 0; i < 256; i++) { j = i; count = j & 1; while (j) count += ((j >>= 1) & 1); bits[i] = count; /* the number of set bits in i */ } /* make checkermask have alternate bits set */ checkermask = 2; /* * Do this 32 times in case we are on a 64 bit machine. The lost * time otherwise is highly insignificant. */ for (i = 0; i < 32; i++) checkermask |= (checkermask << 2); return; } void repaint() /* this fixes the window up whenever it is uncovered */ { int i; XSetPlaneMask(display, gc, AllPlanes); drawbutton(quitbutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT, "quit", 1); drawbutton(pausebutton, 0, 0, BUTTONWIDTH, BUTTONHEIGHT, "pause", 1 - 2 * paused); drawbutton(speedbutton, 0, 0, BOUNDWIDTH, 18, "speed", -1); drawbutton(speedbutton, 1 + ((10 - speed) * (BOUNDWIDTH - 2)) / 11, 1, (BOUNDWIDTH - 2) / 11, 16, "", 2); if (canonical) drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT, "canonical", -1); else drawbutton(canonbutton, 0, 0, font_width * 15, BUTTONHEIGHT, "microcanonical", 1); /* draw the algorithm buttons */ for (i = 0; i < 2; i++) drawbutton(algobutton, 0, i * 18, 68, 18, algotext[i], 1 - 2 * (algorithm == i)); /* draw the heater buttons */ for (i = 0; i < 3; i++) drawbutton(heatwindow, 0, i * HEATERHEIGHT / 3, HEATERWIDTH, HEATERHEIGHT / 3, heatertext[i], 1 - 2 * (i == (heater + 2 * cooler))); /* draw the boundary condition buttons */ for (i = 0; i < 4; i++) drawbutton(boundwindow, 0, i * BOUNDHEIGHT / 4, BOUNDWIDTH, BOUNDHEIGHT / 4, boundtext[i], 1 - 2 * (i == boundary)); /* draw thermometer */ drawbutton(window, TLEFT - 14, TTOP - 14, TWDTH + 28, TBOTTOM - TTOP + 34, "", 2); drawbutton(window, TLEFT - 3, TTOP - 2, TWDTH + 6, TBOTTOM - TTOP + 4, "", 3); XSetForeground(display, gcpen, black); XDrawLine(display, window, gcpen, TLEFT - 3, TTOP - 2, TLEFT + TWDTH + 2, TTOP - 2); XDrawLine(display, window, gcpen, TLEFT - 3, TTOP - 2, TLEFT - 3, TBOTTOM + 2); XDrawLine(display, window, gcpen, TLEFT + TWDTH + 2, TBOTTOM + 2, TLEFT + TWDTH + 2, TTOP - 2); XFillArc(display, window, gcpen, TLEFT - (3 * TWDTH) / 2 - 3, TBOTTOM - TWDTH - 3, 4 * TWDTH + 4, 4 * TWDTH + 4, 0, 64 * 360); XSetForeground(display, gcpen, white); XFillArc(display, window, gcpen, TLEFT - (3 * TWDTH) / 2 - 2, TBOTTOM - TWDTH - 2, 4 * TWDTH, 4 * TWDTH, 0, 64 * 360); XSetForeground(display, gcpen, darkcolor); XFillArc(display, window, gcpen, TLEFT - (3 * TWDTH) / 2, TBOTTOM - TWDTH, 4 * TWDTH - 2, 4 * TWDTH - 2, 0, 64 * 360); XSetForeground(display, gcpen, spinup); XFillArc(display, window, gcpen, TLEFT - TWDTH / 2, TBOTTOM, TWDTH + 2, TWDTH + 2, 10, 64 * 250); XSetForeground(display, gcpen, black); for (i = TBOTTOM - 10; i > TTOP; i -= (THGT / 20)) { XDrawLine(display, window, gcpen, TLEFT - 7, i, TLEFT - 4, i); XDrawLine(display, window, gcpen, TLEFT + TWDTH + 3, i, TLEFT + TWDTH + 6, i); } i = TBOTTOM - THGT * (-0.6 + 1.0 / log(1. + sqrt(2.))); XDrawLine(display, window, gcpen, TLEFT - 15, i, TLEFT + TWDTH + 15, i); XDrawString(display, window, gcpen, TLEFT + TWDTH + 15, i + 4, "T", 1); XDrawString(display, window, gcpen, TLEFT + TWDTH + 23, i + 8, "c", 1); fixthermometer(); /* write various strings */ sprintf(stringbuffer, "%d by %d lattice", ncols, nrows); XDrawString(display, window, gcpen, 150, 92, stringbuffer, strlen(stringbuffer)); XDrawString(display, window, gcpen, 20, 92, "beta=", 5); XDrawString(display, window, gcpen, LEFTOFFSET, 114, "spins:", 6); XDrawString(display, window, gcpen, LEFTOFFSET, 142 + nrows, "changes:", 8); XDrawString(display, window, gcpen, TLEFT - 8, TBOTTOM + 40, "MJC", 3); /* draw border and redraw images */ drawbutton(window, LEFTOFFSET - 4, TOPOFFSET - 4, ncols + 8, nrows + 8, NULL, 4); drawbutton(window, LEFTOFFSET - 4, TOPOFFSET + 28 + nrows - 4, ncols + 8, nrows + 8, NULL, 4); XPutImage(display, window, gc, spinimage, 0, 0, LEFTOFFSET, TOPOFFSET, ncols, nrows); XPutImage(display, window, gc, clusterimage, 0, 0, LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows); /* this might speed things up a little bit */ XSetPlaneMask(display, gc, spinup ^ spindown); return; } void refreshdemons(int step) /* * refresh cluster demons at beta note: Calling this after each * cluster update gives the canonical algorithm. */ /* * To generate a bit set with probability p: write p as a binary * fraction. Compare a random bit string with the digits of this * fraction and find the first place at which they are the same. * Output the bit that occupies that place. This is done here in * parallel on WORDSIZE demons at a time. */ { long morebits, randombits, bitprob, currentbit; double bfactor, prob; long i; int bit, j; bfactor = exp(-beta); for (bit = 0; bit < DEMONBITS; bit++) { bfactor *= bfactor; /* gives exp(-2 beta) for the first * bit, exp(-4 beta) for the next, * ... */ prob = bfactor / (1. + bfactor); /* prob that current demon * bit is 1 */ bitprob = ((1 << 16) * prob); /* convert to 16 bit integer */ for (i = 0; i < ndemons; i += step) { demon[bit][i] = 0; morebits = (-1); for (j = 0; j < 16; j++) { randombits = RND; /* new set of random bits */ currentbit = -((bitprob >> (15 - j)) & 1); /* * set demon to current bit if it equals randombits and has * not been fixed yet */ demon[bit][i] |= (currentbit & randombits & morebits); /* shut off morebits where they are equal */ morebits &= (currentbit ^ randombits); if (0 == morebits) break; /* all bits decided */ } } } return; } void heatbath() /* * refresh local demons at beta */ /* * To generate a bit set with probability p: write p as a binary * fraction. Compare a random bit string with the digits of this * fraction and find the first place at which they are the same. * Output the bit that occupies that place. This is done here in * parallel on WORDSIZE demons at a time. */ { long morebits, randombits, currentbit; int j; demon0 = demon1 = 0; morebits = (-1); for (j = 15; (j >= 0) && (morebits != 0); j--) { randombits = RND; /* new set of random bits */ currentbit = -((bitprob0 >> j) & 1); /* current bit is a word filled with the j'th bit of bitprob0 */ /* * set demon to current bit if it equals randombits and has not * been fixed yet */ demon0 |= (currentbit & randombits & morebits); /* shut off morebits where they are equal */ morebits &= (currentbit ^ randombits); if (0 == morebits) break; /* all bits decided */ } /* now do demon1 */ morebits = (-1); for (j = 15; (j >= 0) && (morebits != 0); j--) { randombits = RND; currentbit = -((bitprob1 >> j) & 1); demon1 |= (currentbit & randombits & morebits); morebits &= (currentbit ^ randombits); if (0 == morebits) break; } return; } void localupdate() /* the local updating routine */ { long *up, *down; /* pointers to above and below rows */ long temp0, temp1, reject, spins, left, right, top, bottom; int row = 0, col, rowloop; checkermask = (~checkermask); for (rowloop = 0; rowloop < nrows; rowloop++) { row = (row + 29); /* jump ahead 29 rows; allows for * energy to move long distances */ while (row >= nrows) row -= nrows; /* get neighbors */ up = field[previousrow]; down = field[nextrow]; leftshift(field[row], work0); rightshift(field[row], work1); /* to keep speed up, only refresh occasionally */ if ((canonical) && ((row & 3) == 3)) heatbath(); for (col = 0; col < hwords; col++) { spins = field[row][col]; /* determine antiparallel neighbors */ right = work0[col] ^ spins; left = work1[col] ^ spins; top = up[col] ^ spins; bottom = down[col] ^ spins; if (0 == row) switch (boundary) { case 1: top ^= (~0); break; case 2: top = (~0) ^ spins; break; case 3: top = spins; break; } else if (nrows - 1 == row) switch (boundary) { case 1: bottom ^= (~0); break; case 2: bottom = (~0) ^ spins; break; case 3: bottom = spins; break; } /* add up top, bottom, and left antiparallel spins */ temp0 = top ^ bottom ^ left; temp1 = (top & (bottom | left)) | (bottom & left); /* add in right */ reject = temp0 & temp1 & right; temp1 ^= (temp0 & right); temp0 ^= right; /* add demon0 */ reject ^= (temp0 & temp1 & demon0); temp1 ^= (temp0 & demon0); temp0 ^= demon0; /* add in demon1 */ reject ^= (temp1 & demon1); temp1 ^= demon1; /* subtract two */ temp1 = ~temp1; reject ^= temp1; /* only change one color sites */ reject |= checkermask; /* make changes */ demon0 = (demon0 & reject) | (temp0 & (~reject)); demon1 = (demon1 & reject) | (temp1 & (~reject)); field[row][col] = (spins ^ (~reject)); /* put changes into cluster image */ cluster[row][col] &= checkermask; cluster[row][col] |= (~reject); energycount += bitcount(demon0 & (~checkermask)); } /* heat or cool */ if (!canonical) if (heater) demon0 |= (RND & RND & RND & RND & RND & RND); else if (cooler) demon1 &= (RND | RND | RND | RND); } /* heat or cool for canonical case */ if (canonical && (heater || cooler)) { beta += (cooler - heater) * .000005 * nrows; bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); fixthermometer(); } /* monitor temperature every 50 "half" iterations */ iteration++; if (iteration == 50) { /* update beta and thermometer */ if (!canonical) { beta = .25 * log((double) (-1. + 50 * volume / (2. * energycount))); bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); fixthermometer(); } iteration = 0; energycount = 0; } /* copy lattice to window */ XPutImage(display, window, gc, spinimage, 0, 0, LEFTOFFSET, TOPOFFSET, ncols, nrows); XPutImage(display, window, gc, clusterimage, 0, 0, LEFTOFFSET, TOPOFFSET + 28 + nrows, ncols, nrows); return; } void fixthermometer() /* fill the thermometer with fluid */ { int temp; sprintf(stringbuffer, "%.3f", beta); XDrawImageString(display, window, gcpen, LEFTOFFSET, 92, stringbuffer, strlen(stringbuffer)); temp = THGT * (-0.6 + 0.5 / beta); if (temp < 0) temp = 0; if (temp > THGT) temp = THGT; XSetForeground(display, gcpen, darkcolor); XFillRectangle(display, window, gcpen, TLEFT, TBOTTOM - temp, TWDTH, temp); XFillRectangle(display, window, gcclear, TLEFT, TTOP, TWDTH, THGT - temp); XSetForeground(display, gcpen, black); return; } void rightshift(source, destination) /* shifts a row of cells to the right by one bit */ long *source, *destination; { int col; long carry = 0; for (col = bottomcol; col != (topcol + nshift); col += nshift) { destination[col] = (~LEFTBIT) & (source[col] >> 1); if (carry) destination[col] |= LEFTBIT; carry = source[col] & 1; } if (carry) destination[bottomcol] |= LEFTBIT; switch (boundary) { case 1: destination[bottomcol] ^= LEFTBIT; break; case 2: destination[bottomcol] |= LEFTBIT; break; case 3: destination[bottomcol] &= (~LEFTBIT); break; } return; } void leftshift(source, destination) /* shifts a row of cells one bit to the left */ long *source, *destination; { int col; long carry = 0; for (col = topcol; col != (bottomcol - nshift); col -= nshift) { destination[col] = (source[col] << 1); if (carry) destination[col] |= 1; carry = source[col] >> (WORDSIZE - 1); } if (carry) destination[topcol] |= 1; switch (boundary) { case 1: destination[topcol] ^= 1; break; case 2: destination[topcol] |= 1; break; case 3: destination[topcol] &= (~1); break; } return; } void setalg() /* take preparatory action if algorithm button hit */ { int row, col; energycount = volumecount = iteration = 0; if (1 == algorithm) { /* reset cluster stuff */ if (boundary > 1) { setboundary(0); } refreshdemons(1); clustergrowth = 0; for (col = 0; col < hwords; col++) for (row = 0; row < nrows; row++) { cluster[row][col] = 0; unchecked[row][col] = 1; } } return; } void setheater(value) /* acts on heater buttons */ int value; { switch (value) { case 0: heater = cooler = 0; break; case 1: heater = 1 - heater; break; case 2: cooler = 1 - cooler; } drawbutton(heatwindow, 0, 0, HEATERWIDTH, HEATERHEIGHT / 3, heatertext[0], 2 * (heater | cooler) - 1); drawbutton(heatwindow, 0, (HEATERHEIGHT) / 3, HEATERWIDTH, HEATERHEIGHT / 3, heatertext[1], 1 - 2 * heater); drawbutton(heatwindow, 0, (HEATERHEIGHT * 2) / 3, HEATERWIDTH, HEATERHEIGHT / 3, heatertext[2], 1 - 2 * cooler); return; } void setboundary(value) /* take action if boundary buttons hit */ int value; { drawbutton(boundwindow, 0, boundary * BOUNDHEIGHT / 4, BOUNDWIDTH, BOUNDHEIGHT / 4, boundtext[boundary], 1); /* don't allow fixed boundaries with cluster algorithm */ if ((algorithm == 1) && (value > 1)) fprintf(stderr, "fixed boundaries not implemented for cluster algorithm\n"); else boundary = value; drawbutton(boundwindow, 0, boundary * BOUNDHEIGHT / 4, BOUNDWIDTH, BOUNDHEIGHT / 4, boundtext[boundary], -1); return; } void openwindow(argc, argv) /* * a lot of this is taken from the basicwin program in the Xlib * Programming Manual */ int argc; char **argv; { char *window_name = "2-d Ising Model"; char *icon_name = "Ising"; Pixmap icon_pixmap; char *display_name = NULL; XColor xcolor, colorcell; Colormap cmap; /* icon */ #define icon_bitmap_width 16 #define icon_bitmap_height 16 static unsigned char icon_bitmap_bits[] = { 0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff}; /* open up the display */ if ((display = XOpenDisplay(display_name)) == NULL) { fprintf(stderr, "%s: cannot connect to X server %s\n", progname, XDisplayName(display_name)); exit(-1); } screen = DefaultScreen(display); depth = DefaultDepth(display, screen); cmap = DefaultColormap(display, screen); /* depth=1; *//* uncomment to test monochrome display */ spindown = darkcolor = black = BlackPixel(display, screen); spinup = lightcolor = white = WhitePixel(display, screen); if (depth > 1) { /* color? This is not the right way * to do it, but ... */ if (XAllocNamedColor(display, cmap, "coral", &colorcell, &xcolor)) darkcolor = colorcell.pixel; if (XAllocNamedColor(display, cmap, "turquoise", &colorcell, &xcolor)) lightcolor = colorcell.pixel; if (XAllocNamedColor(display, cmap, "grey", &colorcell, &xcolor)) spinup = colorcell.pixel; if (XAllocNamedColor(display, cmap, "black", &colorcell, &xcolor)) spindown = colorcell.pixel; } /* make the main window */ windowwidth = ncols + 75; windowheight = 2 * nrows + 164 + 24; window = XCreateSimpleWindow(display, RootWindow(display, screen), 0, 0, windowwidth, windowheight, 4, BlackPixel(display, screen), lightcolor); /* make the icon */ icon_pixmap = XCreateBitmapFromData(display, window, icon_bitmap_bits, icon_bitmap_width, icon_bitmap_height); size_hints.flags = PPosition | PSize | PMinSize; size_hints.min_width = windowwidth; size_hints.min_height = 210; #ifdef X11R3 size_hints.x = 0; size_hints.y = 0; size_hints.width = windowwidth; size_hints.height = windowheight; XSetStandardProperties(display, window, window_name, icon_name, icon_pixmap, argv, argc, &size_hints); #else { XWMHints wm_hints; XClassHint class_hints; XTextProperty windowName, iconName; if (XStringListToTextProperty(&window_name, 1, &windowName) == 0) { fprintf(stderr, "%s: structure allocation for windowName failed.\n" ,progname); exit(-1); } if (XStringListToTextProperty(&icon_name, 1, &iconName) == 0) { fprintf(stderr, "%s: structure allocation for iconName failed.\n" ,progname); exit(-1); } wm_hints.initial_state = NormalState; wm_hints.input = True; wm_hints.icon_pixmap = icon_pixmap; wm_hints.flags = StateHint | IconPixmapHint | InputHint; class_hints.res_name = progname; class_hints.res_class = "Basicwin"; XSetWMProperties(display, window, &windowName, &iconName, argv, argc, &size_hints, &wm_hints, &class_hints); } #endif /* make a default cursor */ XDefineCursor(display,window, XCreateFontCursor(display,XC_hand2)); /* pick the events to look for */ event_mask = ExposureMask | ButtonPressMask | StructureNotifyMask; XSelectInput(display, window, event_mask); /* pick font: 9x15 is supposed to almost always be there */ if ((font = XLoadQueryFont(display, "9x15")) == NULL) { fprintf(stderr, "%s: Cannot open 9x15 font\n", progname); exit(-1); } font_height = font->ascent + font->descent; font_width = font->max_bounds.width; /* * make graphics contexts: gc for black on white, gcpen for * various, defaults to black on lightcolor, gclear does the * reverse, clears thermometer top */ gc = XCreateGC(display, window, 0, NULL); XSetFont(display, gc, font->fid); XSetForeground(display, gc, spinup); XSetBackground(display, gc, spindown); gcpen = XCreateGC(display, window, 0, NULL); XSetFont(display, gcpen, font->fid); XSetForeground(display, gcpen, black); XSetBackground(display, gcpen, lightcolor); gcclear = XCreateGC(display, window, 0, NULL); XSetFont(display, gcclear, font->fid); XSetForeground(display, gcclear, spinup); XSetBackground(display, gcclear, black); /* show the window */ XMapWindow(display, window); return; } void allocarray(array, xsize, ysize, datasize) char ***array; int xsize, ysize, datasize; /* * dynamically allocates a two dimensional (*array)[xize][ysize] with * datasize bytes per element */ { int i; if ((*array) != ((char **) NULL)) /* free previous allocation */ free(*array); /* allocate memory for array of xsize pointers plus data */ *array = (char **) malloc(sizeof(char *) * xsize + datasize * xsize * ysize); if (*array == (char **) NULL) { fprintf(stderr, "%s: memory allocation problems\n", progname); exit(-1); } /* initialize the pointers */ for (i = 0; i < xsize; i++) *((*array) + i) = ((char *) (*array + xsize) + (i * ysize * datasize)); return; } void makebuttons() /* * creates the buttons and allocates arrays; call after any resizing */ { int row, col; XDestroySubwindows(display, window); XSetForeground(display, gcpen, lightcolor); XFillRectangle(display, window, gcpen, 0, 0, windowwidth, windowheight); quitbutton = makebutton(4, 4, BUTTONWIDTH, BUTTONHEIGHT); pausebutton = makebutton(58, 4, BUTTONWIDTH, BUTTONHEIGHT); boundwindow = makebutton(windowwidth - BOUNDWIDTH - 4, 4, BOUNDWIDTH, BOUNDHEIGHT); heatwindow = makebutton(windowwidth - HEATERWIDTH - BOUNDWIDTH - 12, 4, HEATERWIDTH, HEATERHEIGHT); algobutton = makebutton(18, 36, 68, 36); speedbutton = makebutton(windowwidth - BOUNDWIDTH - 4, BOUNDHEIGHT + 26, BOUNDWIDTH, 18); canonbutton = makebutton(LEFTOFFSET + 4 + (ncols - font_width * 15) / 2, windowheight - 26, font_width * 15, BUTTONHEIGHT); /* reset various lattice dimensions */ hwords = ncols / WORDSIZE; ncols = hwords * WORDSIZE; /* truncate ncols to a multiple of * wordsize */ if (0 == (nrows % 29)) nrows--; /* so demon jump of 29 rows works ok */ volume = nrows * ncols; ndemons = 2 * nrows * hwords; /* allocate various arrays */ if (NULL != spinimage) { (*spinimage).data = NULL; XDestroyImage(spinimage); } if (NULL != clusterimage) { (*clusterimage).data = NULL; XDestroyImage(clusterimage); } allocarray(&field, nrows, hwords, sizeof(long)); allocarray(&cluster, nrows, hwords, sizeof(long)); allocarray(&sadx, nrows, hwords, sizeof(long)); allocarray(&sady, nrows, hwords, sizeof(long)); allocarray(&unchecked, nrows, hwords, sizeof(long)); allocarray(&demon, DEMONBITS, ndemons, sizeof(long)); allocarray(&newdemon, DEMONBITS, ndemons, sizeof(long)); if (NULL != work0) free(work0); if (NULL != work1) free(work1); if (NULL != activerow) free(activerow); work0 = (long *) malloc(sizeof(long) * hwords); work1 = (long *) malloc(sizeof(long) * hwords); activerow = (int *) malloc(sizeof(int) * nrows); if ((NULL == work0) || (NULL == work1)) { fprintf(stderr, "%s: memory allocation problems\n", progname); exit(-1); } energycount = volumecount = iteration = 0; /* make image structures */ spinimage = XCreateImage(display, (Visual *) & window, 1, XYBitmap, 0, (char *) *field, ncols, nrows, 32, 0); clusterimage = XCreateImage(display, (Visual *) & window, 1, XYBitmap, 0, (char *) *cluster, ncols, nrows, 32, 0); if ((NULL == spinimage) || (NULL == clusterimage)) { fprintf(stderr, "%s: memory allocation problems\n", progname); exit(-1); } (*spinimage).bitmap_unit = WORDSIZE; (*clusterimage).bitmap_unit = WORDSIZE; /* test for byte order on client */ **field = 1; if (*(char *) (*field)) { (*spinimage).byte_order = LSBFirst; (*clusterimage).byte_order = LSBFirst; /* printf("byte_order=LSBFirst\n"); */ } else { (*spinimage).byte_order = MSBFirst; (*clusterimage).byte_order = MSBFirst; /* printf("byte_order=MSBFirst\n"); */ } /* find out in what order bits are put on screen */ **field = 0xff; if (XGetPixel(spinimage, 0, 0)) { /* bits are backward */ nshift = (-1); topcol = 0; bottomcol = hwords - 1; /* printf("bitmap_bit_order=LSBFirst\n"); */ } else { nshift = 1; topcol = hwords - 1; bottomcol = 0; } /* set initial state partially random and clear cluster array */ for (col = 0; col < hwords; col++) for (row = 0; row < nrows; row++) { field[row][col] = ~(RND * (drand48() > .77)); cluster[row][col] = 0; unchecked[row][col] = 1; } /* equilibrate initial demons at critical temperature */ beta = 0.5 * log(1.0 + sqrt(2.0)); /* critical coupling */ bitprob0 = (1 << 16) * exp(-4 * beta) / (1. + exp(-4 * beta)); bitprob1 = (1 << 16) * exp(-8 * beta) / (1. + exp(-8 * beta)); refreshdemons(1); heatbath(); return; } Window makebutton(xoffset, yoffset, xsize, ysize) int xoffset, yoffset, xsize, ysize; /* * Puts button of specified dimensions on window. Nothing is drawn. */ { Window buttonwindow; long event_mask; Cursor cursor; buttonwindow = XCreateSimpleWindow(display, window, xoffset, yoffset, xsize, ysize, 0, black, lightcolor); event_mask = ButtonPressMask | ExposureMask; /* look for mouse-button * presses */ XSelectInput(display, buttonwindow, event_mask); /* use a hand cursor for buttons */ cursor=XCreateFontCursor(display,XC_hand2); XDefineCursor(display,buttonwindow,cursor); XMapWindow(display, buttonwindow); return buttonwindow; } void drawbutton(buttonwindow, xoffset, yoffset, xsize, ysize, text, state) Window buttonwindow; int xoffset, yoffset, xsize, ysize, state; char *text; /* * Draw a button in buttonwindow of specified dimensions with text * centered. Color of button determined by sign of "state." size of * border determined by magnitude. */ { int textlength, i, j; int cdark, clight, cup, cdown; int cleft, cright, cbutton, ctext; cup = lightcolor; cdown = darkcolor; cdark = black; clight = white; if (state < 0) { cbutton = cdown; ctext = clight; cleft = cdark; cright = clight; } else { cbutton = cup; ctext = cdark; cleft = clight; cright = cdark; } if (1 == depth) clight = cleft = cright = black; j = abs(state); XSetForeground(display, gcpen, cbutton); XFillRectangle(display, buttonwindow, gcpen, xoffset + j, yoffset + j, xsize - 2 * j, ysize - 2 * j); XSetForeground(display, gcpen, cleft); XFillRectangle(display, buttonwindow, gcpen, xoffset, yoffset, xsize, j); XFillRectangle(display, buttonwindow, gcpen, xoffset, yoffset, j, ysize); XSetForeground(display, gcpen, cright); for (i = 0; i < j; i++) { XDrawLine(display, buttonwindow, gcpen, xoffset + i, yoffset + ysize - i - 1, xoffset + xsize - i - 1, yoffset + ysize - i - 1); XDrawLine(display, buttonwindow, gcpen, xoffset + xsize - i - 1, yoffset + i, xoffset + xsize - i - 1, yoffset + ysize - i - 1); } if (NULL != text) { textlength = strlen(text); XSetForeground(display, gcpen, ctext); XDrawString(display, buttonwindow, gcpen, xoffset + (xsize - textlength * font_width) / 2, yoffset + (ysize + font_height) / 2 - 2, text, textlength); } XSetForeground(display, gcpen, darkcolor); return; }