c语言程序ising算法,基于Wolff算法的二维Potts模型蒙特卡罗模拟

8种机械键盘轴体对比

本人程序员，要买一个写代码的键盘，请问红轴和茶轴怎么选？

由于我的本科毕业设计的课题是利用蒙特卡罗模拟的方法求解二维Potts模型的临界温度，所以其中涉及到利用Wolff算法对于Potts模型进行模拟。毕业答辩完成之后将所使用的代码放于此处，便于交流学习和日后巩固。

Wolff 算法

蒙特卡罗模拟有很多的算法，其中比较常用的有Metropolis算法和Wolff算法。由于Metropolis算法在临界点附近效率不高，所以选择Wolff算法来进行模拟。

Wolff算法的主要思想

本段内容来自于统计物理学中的蒙特卡罗方法 1一书。

When the correlation length becomes large close to $T_c$, large regions are often called domains (refers to a group of adjacent spins which are pointing in he same direction). It is quite difficult for the algorithm to flip over one of these large domains because it has to do it spin by spin, and each move has quite a high probability of being rejected. The chances of flipping over a spin in the middle of a domain are particularly low, because it is surrounded by four others pointing in the same direction.

This algorithm was proposed in 1989 by Ulli Wolff.

The basic idea is to look for clusters of similarly oriented spins and then flip them in their entirety all in one go, rather than trying to turn them over spin by painful spin.How to find these clusters of spins which we are going to flip?

Pick a spin at random from the lattice and then look at its neighbours to see if any of them are pointing in the same direction. Then we can look at the neighbours of those neighbours, and so on, iteratively, until we have built up an entire cluster of spins. As we know, at high temperatures, the correlation between spins is weak, and we just need to flip a little number of spins. As we approach $T_c$, we know that the sizes of the clusters become much larger. Therefore, the size of the clusters we flip should increase with decreasing temperature, so we have some probability P of adding a spin, which goes up as the temperature falls.What’s the number of p ?

The precise statement for Potts goes like this: The precise statement for Ising goes like this:Choose a seed spin at random from the lattice.

Look in turn at each of the neighbours of that spin. If they have the same value as the seed spin, add them to the cluster with probability $p=1-e^{-beta J}$.

For each spin that was added in the last step,examine each of its neighbours to find which ones, if any, have the same value and add each of them to the cluster with the same probability P. This step is repeated as many times as necessary until there are no spin left in the cluster whose neighbours have not been considered for inclusion in the cluster.

Choose at random a new value for the spins in the cluster, different from the present value, and set all the spins to that new value.

C语言实现

step() 函数用于完成一次Wolff算法翻转，其中padd, XNN, YNN等变量定义在主函数中。void step()

{

int i;

int sp;

int oldspin, newspin;

int current, nn;

int stack[N];

/* Choose the seed spin for the cluster,

* put it on the stack, and flip it

*/

i = N*genrand_real3();

stack[0] = i;

sp = 1;

oldspin = s[i];

for( ;(newspin = (int)(Q*genrand_real3()))==oldspin; )

{}

s[i] = newspin;

while(sp){

/* Pull a site off the stack */

current = stack[--sp];

/* Check the neighbors */

if((nn=current+XNN)>=N)

nn -= N;

if(s[nn]==oldspin)

if(genrand_real1()

stack[sp++] = nn;

s[nn] = newspin;

}

if((nn=current-XNN)<0)

nn += N;

if(s[nn]==oldspin)

if(genrand_real1()

stack[sp++] = nn;

s[nn] = newspin;

}

if((nn=current+YNN)>=N)

nn -=N;

if(s[nn]==oldspin)

if(genrand_real1()

stack[sp++] = nn;

s[nn] = newspin;

}

if((nn=current-YNN)<0)

nn += N;

if(s[nn]==oldspin)

if(genrand_real1()

stack[sp++] = nn;

s[nn] = newspin;

}

}

}

求解Potts模型的临界温度

并行计算

对于统计模拟，多次重复计算可以降低随机涨落带来的影响。由于本课题利用了大规模集群计算系统，所以需要实现多核并行计算功能，以提高计算效率。具体代码为：int npes,myrank;

MPI_Init(NULL,NULL);

MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

MPI_Comm_size(MPI_COMM_WORLD, &npes);

unsigned long init[4]={0x523, 0x234, 0x345, 0x456}, length=4;

for(int j=0;j<4;j++){

init[j]+=(time(NULL)+myrank);

}

init_by_array(init, length);

for(n=0; n

s[n] = (int)(Q*genrand_real3());

char filename[30]={0},num[5]={0};

strcpy(filename,"file");

sprintf(num, "%d", myrank);

strcat(filename,num);

strcat(filename,".txt");

ofstream outfile(filename);

前四行代码是并行计算所需要的常规设置，中间利用了time()函数和每个核的rank值给格点的自旋方向s[n]进行随机赋值，其中使用到了mt19937ar.h头文件。

后六行代码是希望每个核计算完成后将所算得的数据存在一个名为 filenamemyrank.txt 的文件中，其中myrank是该核在计算开始时被分配到的数字。

求解矩阵本征值

因为求解Potts模型临界温度的过程中需要求解所定义的相似度矩阵的本征值，所以实现求解矩阵本征值的功能。代码如下：MatrixXd Name = MatrixXd::Constant(MC_strps_persite/SpanTime,MC_strps_persite/SpanTime,0.0);

//VectorXcd v;

for(k = 0;k

for(l=0; l

Name(k,l) = Connectmatrix[k][l];

}

SelfAdjointEigenSolver es(Name);

//cout<

//cout<

double eigen_values[MC_strps_persite/SpanTime];

//double eigen_vactors[MC_strps_persite/SpanTime][MC_strps_persite/SpanTime];

for(k=0; k

//v = es.eigenvectors().col(MC_strps_persite/SpanTime-k-1);

eigen_values[k] = es.eigenvalues()(MC_strps_persite/SpanTime-k-1);

//eigen_vactors[k][l] = v(l);

}

//fprintf(fpt_2, "The temperature is %lfn", Temperature);

for(k=0; k

outfile <

outfile <

此处使用了c++中求解本征值的eigen包，由于该包需要使用MatrixXd形数据，所以代码前五行将原始的数组形式的矩阵转换成MatrixXd形的数据。其中注释掉的部分可用于求解特征向量。最后该代码输出了该矩阵前Q*3个本征值。(Q为Potts模型方向数，本征值按从大到小输出)

1: Newman M E J, Barkema G T. Monte Carlo Methods in Statistical Physics chapter 1-4[M]. Oxford University Press: New York, USA, 1999.