Main Program: Monte Carlo Loop

And now the loop. Note the outer loop counts cycles, and the inner loops counts flip attempts in one cycle. The variables i and j are randomly assigned between 0 and $L$-1, identifying a random spin. This random spin is tagged and sent to our dE() function to calculate the change in energy we would expect if that spin were flipped (+1 $\rightarrow$-1 or -1 $\rightarrow$+1). We then calculate the Boltzmann factor in b, and then use the Metropolis criterion to decide whether to accept or reject this spin flip: choosing a random variable on [0,1] and accepting the move if the Boltzmann factor is greater than this number. If it is accepted, we actually peform the flip by multiplying the spin value by -1. Finally, if the current cycle is one in which we collect a sample, we do so by calling the samp() function. The logical expression (c%fSamp)! evaluates to TRUE if the cycle counter c is divisible by fSamp. In the call to samp(), the arguments s and e are passed by reference, signified by the & qualifier. This is necessary because inside samp() we modify both variables. Each of s and e are added to their respective tallies, ssum and esum, and the number of samples taken nSamp is incremented by 1.

  s = 0.0;
  e = 0.0;
  nSamp = 0;
  for (c=0;c<nCycles;c++) {
    /* Make N flip attempts */
    for (a=0;a<N;a++) {
      /* randomly select a spin */
      i=(int)gsl_rng_uniform_int(r,L);
      j=(int)gsl_rng_uniform_int(r,L);
      /* get the "new" energy as the incremental change due
         to flipping spin (i,j) */
      de = dE(M,L,i,j);
      /* compute the Boltzmann factor; recall T is now
         reciprocal temperature */
      b = exp(de*T);
      /* pick a random number between 0 and 1 */
      x = gsl_rng_uniform(r);
      /* accept or reject this flip */
      if (x<b) { /* accept */
	      /* flip it */
	      M[i][j]*=-1;
      }
    }
    /* Sample and accumulate averages */
    if (!(c%fSamp)) {
      samp(M,L,&s,&e);
      fprintf(stdout,"%i %.5le %.5le\n",c,s,e);
      fflush(stdout);
      ssum+=s;
      esum+=e;
      nSamp++;
    }
  }