以下のプログラムは、イジングモデルの数値計算例 1 のプログラムをさらにエネルギーと比熱と磁化を求めるために改良したものである。

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

#define n 100                   //系の最大数
#define J 1.0                   //交換相互作用
#define max 500                 //平均エネルギーを求めるための試行回数

double ising(int seed,double T,double *M);              //イジングモデルの基本的計算を行う関数
double E(double spin[n]);               //エネルギーを計算する関数
int judge(double beta,double oldE,double newE);         //スピンの向きを変えるのかを判断する関数

int main()
{
        int i,seed;
        double T,U,aveU1,aveU2,C,M,m;

        FILE *output1,*output2,*output3;
        output1=fopen("energy.data","w");
        output2=fopen("hinetsu.data","w");
        output3=fopen("zika.data","w");

/*理論的にはTはkT(ボルツマン定数Ｘ温度)に値する*/
        for(T=0.1;T<10.0;T+=0.1) {
                seed=0;         //ここでは乱数の種をただ１ずつ増やすだけである。他の方法でもよい。
                aveU1=aveU2=0;
                m=0.0;
                for(i=0;i<max;i++) {
                        M=0.0;
                        U=ising(seed,T,&M);
                        aveU1+=U;
                        aveU2+=U*U;
                        m+=M;
                        seed++;
                }
                aveU1/=max;             //これはエネルギーの平均値
                aveU2/=max;             //これはエネルギーの２乗の平均値
                m/=max;                 //磁化の平均値
                C=(aveU2-aveU1*aveU1)/(T*T);    //比熱の計算

                fprintf(output1,"%f %f\n",T,aveU1);
                fprintf(output2,"%f %f\n",T,C);
                fprintf(output3,"%f %f\n",T,m);

        }

        return 0;
}

/*乱数の種と温度と磁化の値を入れる変数のポインタを受け取り、
乱数の種と温度の値に応じてイジングモデルのシミュレーションを行う*/
double ising(int seed,double T,double *M)
{
        int i,m;
        double spin[n];
        double beta,oldE,newE,U;

        beta=1.0/T;

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

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

        for(i=0;i<500;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(i=0;i<n;i++) {
                *M+=spin[i];
        }
        *M=*M/n;

        U=E(spin);              //エネルギーの計算

        return U;
}


//スピン相互作用によるエネルギーの計算
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;
}

エネルギーと比熱と磁化についてのグラフをみてみる。横軸は kT (ボルツマン定数 X 温度) 、縦軸はそれぞれの種類の値の大きさである。

種類	グラフ
エネルギー
比熱
磁化

ここで比熱のグラフに注目すると、横軸の 1 付近で急激に値が大きくなっている。このときエネルギーも急激に増え始め、磁化は消え始めている。この現象を相転移といい、系は秩序状態から無秩序状態へ変化しているのである。

比熱が 0 付近あたりでのふるまいは誤差の影響と思われる。エネルギーのグラフからエネルギーはそのあたりでそれほど変化してないので、比熱は 0 付近では、小さな値であると思われるからだ。