24小时热门版块排行榜    

查看: 1169  |  回复: 0

zyj8119

木虫 (著名写手)

[交流] 【转帖】wolff算法的C代码

CODE:
/* Wolff single-cluster algorithm for the two dimensional nearest neighbor
   Ising model                          By Jian-Sheng Wang, February 1995.
   Compile with  

       cc -O Wolff.c -lm

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

#include
#include
#include
                              /* 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 */

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

/*  The main program, running the Monte Carlo loop, collecting data  */

void main()
{
   int i, mc;
   double e = 0;

   p = 1 - exp( - 2/2.269);

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

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

/* 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. */

void monte_carlo_steps(int n)
{
   int i, k;

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

/*  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. */

void flip(int i, int s0)
{
   int j, nn[Z];

   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);
}

/* 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. */

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

   r = i;
   p = 1 - L;
   q = 1;

   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;
   }
}

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

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 */
}

其中drand48() 为0-1之间均匀分布的随机数

回复此楼
好好学习,天天向上。
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

相关版块跳转 我要订阅楼主 zyj8119 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 有没有仍没收到信息的 +6 德尚中行 2026-08-27 7/350 2026-08-30 12:53 by Bxndjjd
[考研] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 11:52 by l0VvVHGBGRLv
[公派出国] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +6 ASdOkHsho7FD 2026-08-28 8/400 2026-08-30 11:32 by l0VvVHGBGRLv
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +6 ASdOkHsho7FD 2026-08-28 9/450 2026-08-30 11:19 by l0VvVHGBGRLv
[考博] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 G6APbkg8SA6w 2026-08-29 4/200 2026-08-30 07:28 by ZPa0EcMwuECS
[硕博家园] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 gy1nBQXYQJqL 2026-08-29 5/250 2026-08-30 06:55 by ZPa0EcMwuECS
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +7 ASdOkHsho7FD 2026-08-28 11/550 2026-08-30 05:40 by ZPa0EcMwuECS
[考博] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 01:34 by ZPa0EcMwuECS
[基金申请] 为什么到现在没收到通知? +4 tannykie 2026-08-29 4/200 2026-08-29 17:14 by lmz0216
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +9 苏东坡二世 2026-08-23 9/450 2026-08-29 14:41 by hustersqt
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 能否退出参与的面上项目解除限项 +23 koalala 2026-08-24 26/1300 2026-08-26 14:29 by 宝贝虫子
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
信息提示
请填处理意见