Exibir código-fonte para Modelo de ising com magnetização constante e algoritmo de Kawasaki

<source lang = "c">
/* beta = inverse temperature
* prob[] = array of acceptance probabilities
* s[] = lattice of spins with helical boundary conditions
* up[] = list of up-pointing spins
* down[] = list of down-pointing spins
* nup = number of up-pointing spins
* ndown = number of down-pointing spins
* L = constant edge length of lattice
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mc.h"
#define LSIZE 16
#define L2    (LSIZE*LSIZE)
//#define TEMP  2.8


void gnuplot_view(void);
void confere_densidade(long int tempo);
void confere_energia(long int tempo);
void confere_energia_densidade(long int tempo);


int s[L2], viz[L2][4] ;
int up[L2], down[L2] ;
int nup,ndown;
int ET;
int MCS;
double prob[17];
double TEMP;
//double beta=1.0/TEMP;

//
void initialize() {

  int i,n1,n2,n3,n4,L; 
  
  // probs
  for (i=0; i<17; i++) {
    prob[i] = exp(-i/(1.0*TEMP));
  }

  // viz
  L = LSIZE;
  for(i=0;i<L2;i++) {
    n1 = ((i/L -1 +L2)%L)*L+i%L;
    n2 = ((i/L)%L)*L+(i+1+L)%L;
    n3 = ((i/L +1 +L2)%L)*L+i%L;
    n4 = ((i/L)%L)*L+(i-1+L)%L;
    
    viz[i][0] = n1;
    viz[i][1] = n2;
    viz[i][2] = n3;
    viz[i][3] = n4;
  }

  // estado -- lembrar que esta versao esta formulada em termos de s= +- 1 (ising)
  nup = 0;
  ndown = 0;
  for(i=0;i<L2;i++) {
    up[i] = 0;
    down[i] = 0;
    if((rand()/(double) RAND_MAX)<0.25) {
      s[i] = 1;
      up[nup] = i;
      nup++;
    }
    else {
      s[i]=-1;
      down[ndown] = i;
      ndown++;
    }
    
  }

  ET = 0;
  for(i = 0; i < L2; i++){
	ET += (s[viz[i][0]] + s[viz[i][1]] + s[viz[i][2]] + s[viz[i][3]])*s[i];
  }

  ET = -ET/2;

  
  
}

//
void sweep(long int tempo) {
  
  int i;
  int delta;
  int iup,idown;
  int xup,xdown;
  int upnn1,upnn2,upnn3,upnn4;
  int downnn1,downnn2,downnn3,downnn4;
  int term1,term2;
  double fprob,frand,frand1,frand2;
  
  for (i=0; i<L2; i++) {

    /* Choose a pair of spins */
    iup = rand()%nup;
    idown = rand()%ndown; 
    //
    xup = up[iup];
    xdown = down[idown];
    
    upnn1 = viz[xup][0];
    upnn2 = viz[xup][1];
    upnn3 = viz[xup][2];
    upnn4 = viz[xup][3];
    downnn1 = viz[xdown][0];
    downnn2 = viz[xdown][1];
    downnn3 = viz[xdown][2];
    downnn4 = viz[xdown][3];
    
    /* Calculate the local energies before swapping */
    term1 = s[upnn1]+s[upnn2]+s[upnn3]+s[upnn4];
    term2 = -s[downnn1]-s[downnn2]-s[downnn3]-s[downnn4];

    /* Swap the spins over */
    s[xup] = -1;
    s[xdown] = +1;

    /* Calculate the difference in the local energies after */
    // original
    term1 -= -s[upnn1]-s[upnn2]-s[upnn3]-s[upnn4] ;
    // modificado por clareza
    //term1 += s[upnn1]+s[upnn2]+s[upnn3]+s[upnn4];
    // original
    term2 -= s[downnn1]+s[downnn2]+s[downnn3]+s[downnn4];
    // modificado por clareza
    //term2 += -s[downnn1]-s[downnn2]-s[downnn3]-s[downnn4];

    /* Calculate total change in energy */
    delta = term1 + term2;

    /* Accept or reject the move */
    fprob = prob[delta];
    frand = rand()/(double) RAND_MAX;
    
    
    if (delta<=0) {
      up[iup] = xdown;
      down[idown] = xup;
      ET = ET + delta;	// MEDIDA DE ENERGIA
    }
    else if (frand < fprob) {
      up[iup] = xdown;
      down[idown] = xup;
      ET = ET + delta;	// MEDIDA DE ENERGIA
    }
    else {
      s[xup] = +1;
      s[xdown] = -1;
    }

  
}


void sweep_local(void){

	int i;
	int Ei, Ef;
	int delta;
	int iup,idown;
	int xup,xdown;
	int upnn1,upnn2,upnn3,upnn4;
	int downnn1,downnn2,downnn3,downnn4;
	int term1,term2;

	for (i=0; i<L2; i++) {

		/* Choose a pair of spins */
		iup = L2*(rand()/(double) RAND_MAX);
		idown = 4*(rand()/(double) RAND_MAX);
		xup = iup;
		xdown = viz[xup][idown];

		if(s[xup] != s[xdown]){
		
		  upnn1 = viz[xup][0];
		  upnn2 = viz[xup][1];
		  upnn3 = viz[xup][2];
		  upnn4 = viz[xup][3];
		  downnn1 = viz[xdown][0];
		  downnn2 = viz[xdown][1];
		  downnn3 = viz[xdown][2];
		  downnn4 = viz[xdown][3];
		  
		  /* Calculate the local energies before swapping */
		  term1 = s[upnn1]+s[upnn2]+s[upnn3]+s[upnn4];
		  term2 = s[downnn1] + s[downnn2] + s[downnn3] + s[downnn4];
      Ei = -s[xup]*term1 - s[xdown]*term2;
		  
      term1 = term1 + s[xup] - s[xdown];
      term2 = term2 + s[xdown] - s[xup];
		
      
      Ef = -s[xdown]*term1 - s[xup]*term2;

		  /* Calculate total change in energy */
      delta = Ef - Ei;

		  /* Accept or reject the move */
		  if (delta<=0) {
        		s[xup] = -s[xup];
        		s[xdown] = -s[xdown];
			    ET += delta;	// MEDIDA DE ENERGIA
		  }
		  else if ((rand()/(double) RAND_MAX) <prob [delta]) {
  		  		s[xup] = -s[xup];
      			s[xdown] = -s[xdown];
			    ET += delta;	// MEDIDA DE ENERGIA
 		}
		else {

  		

		}

	}

  }

}



//////////////////////////////////////////////////////////	 
int main (int argc, char *argv[]) {

  int steps;
  int SEMENTE;
  int TMAX, VIEW, TRAN;
  FILE *arq;
  char nome_do_arquivo[300];
  double temperatura;

#ifdef DEBUG
  temperatura = 2.8;
  TEMP = temperatura;
#else
  temperatura = atof(argv[1]);
  TEMP = temperatura;
#endif

  TMAX=1E6;
  TRAN=1E5;
  VIEW=1;

  /* seed = start_randomic(); */
#ifdef DEBUG
  SEMENTE = 24242;
  srand(SEMENTE);
#else
  SEMENTE = atoi(argv[2]);
  srand(SEMENTE);
#endif

  initialize();
  confere_energia_densidade(0); 
  
  //TRANSIENTE
  for(steps=0; steps<TRAN; steps++){
#ifndef LOCAL
    sweep(steps);
#else
    sweep_local();
#endif
  }

#ifdef DEBUG
  sprintf(nome_do_arquivo, "series_backup.dat");
#else
  sprintf(nome_do_arquivo, "series_%.3lf_%d.dat", TEMP, SEMENTE);
#endif

  arq = fopen(nome_do_arquivo, "w");
  fprintf(arq, "#mcs ET S%d L%d\n", SEMENTE, LSIZE);
  for(steps=0;steps<TMAX;steps++) {
#ifndef LOCAL
    sweep(steps);
    fprintf(arq, "%d %d\n", steps, ET);
#else
    sweep_local();
#endif
#ifdef GNU
    MCS = steps;
    if(steps > 0) {
      gnuplot_view();
      confere_densidade(steps); 
      confere_energia(steps); 
    }
#endif
    //confere_densidade(steps); 
    //confere_energia(steps);
    confere_energia_densidade(steps); 
  

  }

  return 0;
}


/****************************************************************************
 * 		    	       GNUPLOT VIEW		 	 	    *
 ***************************************************************************/
void gnuplot_view(void) {

  int ii,jj,L;

  L=LSIZE;
  
  //printf("set xrange [0:%d]\nset yrange [0:%d]\n", 2*(L-1), (L-1));
  //printf("set cbrange [0:5]\n");
  printf("unset xtics\nunset ytics\n");
  //printf("set palette defined (0 'black', 1 'dark-violet')\n");
  //printf("unset colorbox\n");
  printf("set size square\n");
  printf("set title \"MCS = %d; ET = %d\"\n", MCS, ET);
  //printf("set size ratio 0.5\n");
  printf("set key off\n");
  printf("plot \"-\" matrix w image\n");
  
  for(jj=0;jj<L;jj++) {
    for(ii=0;ii<L;ii++) {
      printf("%d ", s[ii + jj*L]);
      //	printf("%d ", neigh[ii + jj*L][4]);
    }
    /*printf("%d %d %d ",5,5,5);
    for(ii=0;ii<L;ii++)	{
      //printf("%d ", site[ii + jj*L]);
      printf("%d ", viz[ii + jj*L][4]);
    }
    */
    printf("\n");
    
  }
  
  //printf("\ne\n pause 1.0\n");
  printf("\ne\n pause 0.1\n");
  
  return;
}



/****************************************************************************
 * 		    	       Confere densidade         	 	    *
 ***************************************************************************/
void confere_densidade(long int tempo) {

  int i,dens_up,dens_down;

  dens_up = 0;
  dens_down = 0;
  for(i=0;i<L2;i++) {
    if(s[i]==1) {
      dens_up++;
    }
    else if(s[i]==-1) {
      dens_down++;

    }
  }

  if((dens_up!=nup)||(dens_down!=ndown)||(dens_up+dens_down!=L2)) {
    fprintf(stderr,"ERRO!!! no numero de particulas\n");
    fprintf(stderr,"        tempo:     %ld\n",tempo);
    fprintf(stderr,"        dens_up:   %d  nup:   %d\n",dens_up,nup);
    fprintf(stderr,"        dens_down: %d  ndown: %d\n",dens_down,ndown);
    fprintf(stderr,"        dens:      %d  n:     %d\n",dens_up+dens_down,nup+ndown);
    exit(-1);
  }
  return;
}


/****************************************************************************
 * 		    	       Confere energia           	 	    *
 ***************************************************************************/
void confere_energia(long int tempo) {

  int i,ene;

  ene = 0;
  for(i = 0; i < L2; i++){
    ene += (s[viz[i][0]] + s[viz[i][1]] + s[viz[i][2]] + s[viz[i][3]])*s[i];
  }
  ene = -ene/2;

  if(ene!=ET) {
    fprintf(stderr,"ERRO!!! na energia\n");
    fprintf(stderr,"        tempo: %ld   ene: %d  ET: %d\n\n",tempo,ene,ET);
    exit(-1);
  }
  return;
}



/****************************************************************************
 * 		      Confere energia e densidade           	 	    *
 ***************************************************************************/
void confere_energia_densidade(long int tempo) {

  int i,ene,dens_up,dens_down;

  dens_up = 0;
  dens_down = 0;
  for(i=0;i<L2;i++) {
    if(s[i]==1) {
      dens_up++;
    }
    else if(s[i]==-1) {
      dens_down++;

    }
  }

  ene = 0;
  for(i = 0; i < L2; i++){
    ene += (s[viz[i][0]] + s[viz[i][1]] + s[viz[i][2]] + s[viz[i][3]])*s[i];
  }
  ene = -ene/2;


  if((dens_up!=nup)||(dens_down!=ndown)||(dens_up+dens_down!=L2)||(ene!=ET)) {
    fprintf(stderr,"ERRO!!! \n");
    fprintf(stderr,"        tempo:     %ld\n",tempo);
    fprintf(stderr,"        dens_up:   %d  nup:   %d\n",dens_up,nup);
    fprintf(stderr,"        dens_down: %d  ndown: %d\n",dens_down,ndown);
    fprintf(stderr,"        dens:      %d  n:     %d\n",dens_up+dens_down,nup+ndown);
    fprintf(stderr,"        ene:       %d  ET:    %d\n\n",ene,ET);
    exit(-1);
  }

  
  return;
}
</source>