ChinaNet

标题: [转帖]wolff算法的C代码 [打印本页]

作者: Doublen 时间: 2015-4-8 10:22
标题: [转帖]wolff算法的C代码

[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

欢迎光临 ChinaNet (http://www.nuclear.net.cn/)