[转帖]wolff算法的C代码

Doublen

[size=13.63636302947998px]/* Wolff single-cluster algorithm for the two dimensional nearest neighbor
   Ising model                          By Jian-Sheng Wang, February 1995.
   Compile with

[size=13.63636302947998px]       cc -O Wolff.c -lm

   -O for optimization, and -lm to link the math library.               
*/

[size=13.63636302947998px]#include <stdio.h>
#include <stdlib.h>
#include <math.h>
                              /* macro definitions */
#define  L  16                /* lattice linear size */
#define  N  (L*L)             /* total number of spins */
#define  Z  4                 /* coordination number = 2*d */
#define  MCTOT 500            /* total Monte Carlo steps */
#define  MCDIS 200            /* steps discarded in the beginning */
                              /* global variables */
int s[N];                     /* spin +1 or -1 */
double  p = 1.0;              /* percolation probability */

[size=13.63636302947998px]         /* funcition prototypes */
void neighbor(int i, int nn[ ]);
void flip(int i, int s0);
void monte_carlo_steps(int n);
void energy(double *);

[size=13.63636302947998px]/*  The main program, running the Monte Carlo loop, collecting data  */

[size=13.63636302947998px]void main()
{
   int i, mc;
   double e = 0;

[size=13.63636302947998px]   p = 1 - exp( - 2/2.269);

   for (i = 0; i < N; ++i)     /* initialize, all spin up */
      s = 1;

[size=13.63636302947998px]   for(mc = 0; mc < MCTOT; ++ mc) {
      monte_carlo_steps(5);
      if( mc >= MCDIS)
         energy(&e);
   }
   printf("<e> =  %f\n", e/(MCTOT-MCDIS)/N);
}

[size=13.63636302947998px]/* This function monte_carlo_steps performs n cluster flips, eqivalent to
few Monte Carlo steps in standard single-spin-flip algorithms.
It picks a seed site at random and calls the flip function to generate
one cluster.  N is total number of spin.  They are macro definitions. */

[size=13.63636302947998px]void monte_carlo_steps(int n)
{
   int i, k;

[size=13.63636302947998px]   for(k = 0; k < n; ++k) {
      i = drand48() * (double) N;
      flip(i, s);
   }
}

[size=13.63636302947998px]/*  Perform a Wolff single cluster flip. s[], p, and Z are passed globally.
The first argument i of flip function is the site to be flipped, the
second argument is the spin of the cluster before flipping. */

[size=13.63636302947998px]void flip(int i, int s0)
{
   int j, nn[Z];

[size=13.63636302947998px]   s = - s0;                    /* flip the spin immediately */
   neighbor(i, nn);                /* find nearest neighbor of i */
   for(j = 0; j < Z; ++j)          /* flip the neighbor if ...  */
      if(s0 == s[nn[j]] && drand48() < p)
         flip(nn[j], s0);
}

[size=13.63636302947998px]/* Neighbor returns in the array nn[ ] the neighbor sites of i.  The sites
are labelled sequentially, starting from 0.  It works for any hypercubic
lattice.  Z (=2*D) is the coordination number, passed as a macro defintion.
L is linear size, also passed as a macro definition. */

[size=13.63636302947998px]void neighbor(int i, int nn[ ])
{
   int j, r, p, q;

[size=13.63636302947998px]   r = i;
   p = 1 - L;
   q = 1;

[size=13.63636302947998px]   for(j = 0; j < Z; j += 2) {
      nn[j] = (r + 1) % L == 0 ? i + p : i + q;
      nn[j+1]     = r % L == 0 ? i - p : i - q;
      r = r/L;
      p *= L;
      q *= L;
   }
}

[size=13.63636302947998px]/* This function calculate the energy of the configuration s[],
   in fact, it is the negative of the energy.  */

[size=13.63636302947998px]void energy(double *e)
{
   int i, j, ie = 0;
   int nn[Z];
   
   for(i = 0; i < N; ++i) {
      neighbor(i, nn);             /* find the neighbor of center site */
      for(j = 0; j < Z; j += 2)    /* look at positive direction only */
         ie += s*s[nn[j]];      /* add nearest neighbor interaction */
   }
   *e += ie;                       /* accumulant energy */
}

[size=13.63636302947998px]其中drand48() 为0－1之间均匀分布的随机数

[size=13.63636302947998px]本帖转自52MC