１次元イジングモデル。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define n 500                   //系の最大数
#define tmax 1000
#define J 1.0                   //交換相互作用
#define kT 10.0         //ボルルマン定数×温度
#define seed 1127               //乱数の発生の仕方を変える種

double E(double spin[n]);               //エネルギーを計算する関数
int judge(double beta,double oldE,double newE);         //スピンの向きを変えるのかを判断する関数

int main()
{
        int i,j,m;
        double spin[n];
        double beta,oldE,newE;

        beta=1.0/(kT);

        srand(seed);            //乱数の発生の仕方を変える

        FILE *output1,*output2;
        output1=fopen("spinup.data","w");
        output2=fopen("spindown.data","w");

        for(i=0;i<n;i++) {
                spin[i]=1.0;            //スピン初期状態の決定
        }

        for(i=0;i<tmax;i++) {
                oldE=E(spin);           //前のエネルギー値
                m=(int)1.0/(RAND_MAX+1.0)*rand()*n;
                spin[m]*=-1.0;          //スピンの向きを変える
                newE=E(spin);           //新しいエネルギー

                if(judge(beta,oldE,newE)) spin[m]*=-1.0;                //１が帰ってきたらスピンの方向を戻す

                for(j=0;j<n;j++) {
                        if(spin[j]==1.0) {
                                fprintf(output1,"%d %d\n",i,j);
                        }
                        if(spin[j]==-1.0) {
                                fprintf(output2,"%d %d\n",i,j);
                        }
                }
        }

        fclose(output1);
        fclose(output2);

        return 0;
}

//スピン相互作用によるエネルギーの計算
double E(double spin[n])
{
        int i;
        double E;

        E=0.0;

        for(i=0;i<(n-1);i++) {
                E+=-J*spin[i]*spin[i+1];
        }

        return E;
}

//新しい状態が実現するなら０、しないなら１を返す
int judge(double beta,double oldE,double newE)
{
                if((newE>oldE) && (exp(-beta*(newE-oldE)) <= 1.0/(RAND_MAX+1.0)*rand())) return 1;
                else return 0;
}

横軸が時間経過、縦軸が１次元の系の位置を表す。赤色がスピン上向き、緑色がスピン下向きを表す。 様々な結果が得られるので初期値や条件を変えて数値計算してみるとよい。