Я моделирую модель Поттса в двумерной решетке, используя c

Question

Я просто хочу подтвердить, не будет ли каких-либо колебаний в спиновых состояниях при низкой температуре = 0,1, 0,2 и т. Д. Для модели Поттса с q = 2 состояниями.

Достигнет ли он состояния со всеми одинаковыми спинами? Не покажут ли какие-либо флуктуации спиновых состояний при этих температурах во время моделирования?

Если вы хотите скомпилировать код, добавьте -lgsl -lgslcblas во время компиляции.

#include<stdlib.h>
#include<math.h>
#include<stdio.h>
#include<gsl/gsl_rng.h>
#include<string.h>
#include<time.h>
#include"parameters.h"
enum {N=PAR_NO};
int potts[N][N];
double T=PAR_TEMP;
unsigned long int Seed=time(NULL);
gsl_rng *r;
int q=PAR_STATE;
int R=PAR_RSTATE;
void initialise()
{
 for(int i=0;i<N;i++)
  { 
    for(int j=0;j<N;j++)
      {
         potts[i][j]=gsl_rng_uniform_int(r,q+R)+1;
       }
  }
}
double neighbours(int dir,int i,int j)
{
  double neighbour=0;
  if(dir==1)
   {
    if(i==N-1)
      i=0;
    else
     i=i+1;

   }
  if(dir==2)
   {
     if(i==0)
      i=N-1;
     else
      i=i-1;

    }
  if(dir==3)
   { 
     if(j==N-1)
       j=0;
     else
       j=j+1;

   }
  if(dir==4)
   {
       if(j==0)
         j=N-1;
       else
         j=j-1;


    }
 return(potts[i][j]);
}
int codelta(int s,int s1)
{
  if(s==s1)
    return(1);
  else
    return(-1);
}
int delta(int s0,int s1)
{ 
  if(s0==s1)
    return(1);
  else
    return(0);
}
void metropolis(double expo[10])
{
  int a,b,c,dele;
  double ie,fe,s[5];
  for(int i=0;i<N;i++)
  {
    for(int j=0;j<N;j++)
     {
      a=gsl_rng_uniform_int(r,N);
      b=gsl_rng_uniform_int(r,N);
      c=gsl_rng_uniform_int(r,q+R)+1;
      s[0]=potts[a][b];
      s[1]=neighbours(1,a,b);
      s[2]=neighbours(2,a,b);
      s[3]=neighbours(3,a,b);  
      s[4]=neighbours(4,a,b);
      ie=-(delta(s[0],s[1])+delta(s[0],s[2])+delta(s[0],s[3])+delta(s[0],s[4]));
      s[0]=c;
      fe=-(delta(s[0],s[1])+delta(s[0],s[2])+delta(s[0],s[3])+delta(s[0],s[4]));
      dele=fe-ie;
      if(dele>0)
      {
        double x=gsl_rng_uniform(r);
        double P=expo[dele];
        if(x<=P)
         {
           potts[a][b]=s[0];
         }

      }
      else
        potts[a][b]=s[0];
     }

  }

}

double avg_mag()
{ double m1,M[PAR_RSTATE+PAR_STATE],mag,m2,m;
   for(int t=1;t<=q+R;t++)
     {
      for(int i=0;i<N;i++)
       { 
         for(int j=0;j<N;j++)
          {
            if(potts[i][j]==t)
              m1+=1;
            else
              m2+=1;



          }
        }
      M[t]=fabs((m1-m2)/N/N);

       m1=m2=0;
       mag+=M[t];    
     }
     //m=fabs(m1-m2);
   //printf("%lf \n",mag/(q+R));
   return(mag/(q+R));
}

double energy()
{
   double eright,eleft,e;
  for(int i=0;i<N;i++)
   {
    for(int j=0;j<N;j++)
     {
      eright=delta(potts[i][j],potts[i][(j+1)%N]);
      eleft=delta(potts[i][j],potts[(i+1)%N][j]);
      e-=eright+eleft;
     }
   }
  //printf("%lf \n",e);
 return(e);
}
void printmag_ener(double M,double E,int trial,double T,int step)
 {
   char basenm[256];
   sprintf(basenm,"pottsmag_q_%d_r_%d_trial_%d_T_%lf.txt",q,R,trial,T);
   FILE *fp;
   fp=fopen(basenm,"a");
   fprintf(fp,"%d %lf %lf \n",step,E,M);
   fclose(fp);
  }
void print(int step,int trial,double T)
 {
   char basenm[256];
   sprintf(basenm,"potts_q_%d_r_%d_trial_%d_T_%lf_step_%d.txt",q,R,trial,T,step);
   FILE *fp;
   fp=fopen(basenm,"w");
   for(int i=0;i<N;i++)
   { for(int j=0;j<N;j++)
      fprintf(fp,"%d %d %d \n",i,j,potts[i][j]);
   }
   fclose(fp);
  }
int main()
{
 const gsl_rng_type * Typ;
 Typ=gsl_rng_mt19937;
 int trial=PAR_TRIAL;
 r=gsl_rng_alloc(Typ);
 gsl_rng_set(r,Seed);
 double expo[10],E,M;

    for(int i=0;i<10;i++)
      {expo[i]=exp(-i/T);}
     initialise();

    for(int step=0;step<PAR_MCSTEPS;step=step+10)
       {
           printf("%d \n",step);
           M=avg_mag();
           E=energy();
           print(step,trial,T);
           printmag_ener(M,E/N/N,trial,T,step);
          for(int j=0;j<10;j++)
          {
             metropolis(expo);

          }

       }


  return(0);
}

Я ожидаю, что намагниченность будет 1 в идеальном случае, но при симуляции я думаю, что она будет 0,99 или 0,998 или 1 и т. Д. Я имею в виду, разве это не изменится хотя бы немного?