Advertisement
Guest User

Untitled

a guest
Jun 30th, 2012
702
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #define ANISOTROPIC
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #include <stdlib.h>
  6. #include <string.h>
  7.  
  8. #define LATTICE_X   (30)
  9. #define LATTICE_Y   (30)
  10. #define LATTICE_Z   (30)
  11.  
  12. #define MIN(a,b)    (((a)>(b))?(b):(a))
  13.  
  14. struct configuration_t
  15. {
  16.     double phase[LATTICE_X][LATTICE_Y][LATTICE_Z];
  17. };
  18.  
  19. void clear_configuration(struct configuration_t *cfgt)
  20. {
  21.     int i,j,k;
  22.    
  23.     for(i=0;i<LATTICE_X;i++)
  24.         for(j=0;j<LATTICE_Y;j++)
  25.             for(k=0;k<LATTICE_Z;k++)
  26.                 cfgt->phase[i][j][k]=0.0f;
  27. }
  28.  
  29. double gen_random_number(void)
  30. {
  31.     return drand48();
  32. }
  33.  
  34. int gen_random_int(void)
  35. {
  36.     return lrand48();
  37. }
  38.  
  39. #define INIT_TO_ZERO    (0x1)
  40.  
  41. struct configuration_t *init_configuration(short flags)
  42. {
  43.     struct configuration_t *ret=malloc(sizeof(struct configuration_t));
  44.  
  45.     if(!ret)
  46.         return NULL;
  47.    
  48.     if((flags&INIT_TO_ZERO)!=0)
  49.         clear_configuration(ret);
  50.    
  51.     return ret;
  52. }
  53.  
  54. void fini_configuration(struct configuration_t *cfgt)
  55. {
  56.     if(cfgt)
  57.         free(cfgt);
  58. }
  59.  
  60. double lsite_energy(struct configuration_t *cfgt,int x,int y,int z)
  61. {
  62.     double theta1,theta2,xy=0.0f;
  63.  
  64.     theta1=cfgt->phase[x][y][z];
  65.     theta2=cfgt->phase[(x+1)%LATTICE_X][y][z];
  66.  
  67. #ifdef  ANISOTROPIC
  68.     xy-=2.0f*cos(theta1-theta2);
  69. #else
  70.     xy-=cos(theta1-theta2);
  71. #endif
  72.  
  73.     theta1=cfgt->phase[x][y][z];
  74.     theta2=cfgt->phase[x][(y+1)%LATTICE_Y][z];
  75.  
  76.     xy-=cos(theta1-theta2);
  77.  
  78.     theta1=cfgt->phase[x][y][z];
  79.     theta2=cfgt->phase[x][y][(z+1)%LATTICE_Z];
  80.  
  81.     xy-=cos(theta1-theta2);
  82.  
  83.     return xy;
  84. }
  85.  
  86. double configuration_energy(struct configuration_t *cfgt)
  87. {
  88.     int i,j,k;
  89.     double ret=0.0f;
  90.    
  91.     for(i=0;i<LATTICE_X;i++)
  92.         for(j=0;j<LATTICE_Y;j++)
  93.                 for(k=0;k<LATTICE_Z;k++)
  94.                     ret+=lsite_energy(cfgt,i,j,k);
  95.  
  96.     return ret;
  97. }
  98.  
  99. void update_lsite(struct configuration_t *cfgt,int x,int y,int z)
  100. {
  101.    
  102. #define AMPLITUDE   (2*M_PI)
  103.  
  104.     cfgt->phase[x][y][z]+=(gen_random_number()*AMPLITUDE);
  105. }
  106.  
  107. void push_lsite(struct configuration_t *cfgt,int x,int y,int z,double *s)
  108. {
  109.     s[0]=cfgt->phase[x][y][z];
  110. }
  111.  
  112. void pop_lsite(struct configuration_t *cfgt,int x,int y,int z,double *s)
  113. {
  114.     cfgt->phase[x][y][z]=s[0];
  115. }
  116.  
  117. int update_configuration(struct configuration_t *cfgt,double beta)
  118. {
  119.     int i,j,k,x,y,z;
  120.     double deltae,boltzmann_factor,p;
  121.     double f[4];
  122.  
  123.     i=gen_random_int()%LATTICE_X;
  124.     j=gen_random_int()%LATTICE_Y;
  125.     k=gen_random_int()%LATTICE_Z;
  126.  
  127.     deltae=0.0f;
  128.    
  129.     for(x=i-1;x<=(i+1);x++)
  130.         for(y=j-1;y<=(j+1);y++)
  131.             for(z=k-1;z<=(k+1);z++)
  132.                 deltae-=lsite_energy(cfgt,(x+LATTICE_X)%LATTICE_X,(y+LATTICE_Y)%LATTICE_Y,(z+LATTICE_Z)%LATTICE_Z);
  133.  
  134.     push_lsite(cfgt,i,j,k,f);
  135.     update_lsite(cfgt,i,j,k);
  136.  
  137.     for(x=i-1;x<=(i+1);x++)
  138.         for(y=j-1;y<=(j+1);y++)
  139.             for(z=k-1;z<=(k+1);z++)
  140.                 deltae+=lsite_energy(cfgt,(x+LATTICE_X)%LATTICE_X,(y+LATTICE_Y)%LATTICE_Y,(z+LATTICE_Z)%LATTICE_Z);
  141.  
  142.     boltzmann_factor=exp(-beta*deltae);
  143.     p=MIN(1,boltzmann_factor);
  144.  
  145.     if(p<gen_random_number())
  146.     {
  147.         pop_lsite(cfgt,i,j,k,f);
  148.         return 1;
  149.     }
  150.    
  151.     return 2;
  152. }
  153.  
  154. double magnetization(struct configuration_t *cfgt)
  155. {
  156.     int i,j,k;
  157.     double a,b;
  158.    
  159.     a=b=0.0f;
  160.     for(i=0;i<LATTICE_X;i++)
  161.     {
  162.         for(j=0;j<LATTICE_Y;j++)
  163.         {
  164.                 for(k=0;k<LATTICE_Z;k++)
  165.                 {
  166.                     a+=cos(cfgt->phase[i][j][k]);
  167.                     b+=sin(cfgt->phase[i][j][k]);
  168.                 }
  169.         }
  170.     }
  171.  
  172.     return sqrt(pow(a,2.0f)+pow(b,2.0f));
  173. }
  174.  
  175. double magnetization_per_spin(struct configuration_t *cfgt)
  176. {
  177.     double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;
  178.    
  179.     return magnetization(cfgt)/spins;
  180. }
  181.  
  182. double energy_per_spin(struct configuration_t *cfgt)
  183. {
  184.     double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;
  185.    
  186.     return configuration_energy(cfgt)/spins;
  187. }
  188.  
  189. double helicity_modulus(struct configuration_t *cfgt,double T)
  190. {
  191.     int i,j,k;
  192.     double f,g;
  193.     double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;
  194.  
  195.     f=configuration_energy(cfgt);
  196.    
  197.     g=0.0f;
  198.     for(i=0;i<LATTICE_X;i++)
  199.         for(j=0;j<LATTICE_Y;j++)
  200.             for(k=0;k<LATTICE_Z;k++)
  201.                 g+=sin(cfgt->phase[i][j][k]-cfgt->phase[(i+1)%LATTICE_X][j][k]);
  202.  
  203.     return (-0.5f*f-(1.0f/T)*pow(g,2.0f))/spins;
  204. }
  205.  
  206. int main(void)
  207. {
  208.     int c,d;
  209.    
  210.     struct configuration_t *cfgt=init_configuration(INIT_TO_ZERO);
  211.     double T,beta,mpsmean,epsmean,hmmean;
  212.  
  213.     int thermalization=100000;
  214.     int iterations=150000;
  215.     int valid;
  216.     int runs;
  217.  
  218.     for(T=0.0f;T<=50.0f;T+=0.25f)
  219.     {
  220.         epsmean=mpsmean=hmmean=0.0f;
  221.         beta=1.0/T;
  222.         valid=0;
  223.         runs=1;
  224.  
  225.         for(d=0;d<=runs;d++)
  226.         {
  227.             clear_configuration(cfgt);
  228.  
  229.             for(c=0;c<=iterations;c++)
  230.             {
  231.                 double mps,eps,hm;
  232.                 int result=update_configuration(cfgt,beta);
  233.  
  234.                 if((c>thermalization)&&(result==2))
  235.                 {
  236.                     mps=eps=0.0f;
  237.                     mps=magnetization_per_spin(cfgt);
  238.                     eps=energy_per_spin(cfgt);
  239.                     hm=helicity_modulus(cfgt,T);
  240.  
  241.                     mpsmean+=mps;
  242.                     epsmean+=eps;
  243.                     hmmean+=hm;
  244.                     valid++;
  245.                 }
  246.             }
  247.         }
  248.  
  249.         epsmean/=((double)(valid));
  250.         mpsmean/=((double)(valid));
  251.         hmmean/=((double)(valid));
  252.  
  253.         printf("%f %d %f %f %f\n",T,valid,mpsmean,epsmean,hmmean);
  254.         fflush(stdout);
  255.     }
  256.  
  257.     fini_configuration(cfgt);
  258.  
  259.     return 0;
  260. }
Advertisement
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement