Biblioteca com funções em C para simulações de Monte Carlo
Ir para navegação
Ir para pesquisar
/********************************************************************
*** Monte Carlo Routines ***
*** V0.11 18/05/2000 ***
*** ***
*** ***
********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <string.h>
#include <time.h>
#define MC_VERSION "0.11"
#define FNORM (2.3283064365e-10)
#define RANDOM ( (ira[ip++]=ira[ip1++]+ira[ip2++]) ^ira[ip3++] )
#define FRANDOM (FNORM * RANDOM)
/********************************************************************
*** Variable Declarations ***
********************************************************************/
unsigned long randomize, seed; /* Random number */
static clock_t begin_cpu_time; /* Time variables */
static time_t begin_real_time;
unsigned zseed, ira[256];
unsigned char ip,ip1,ip2,ip3;
/********************************************************************
* Random Number Generator by Parisi & Rapuano *
* Last Modified: 18/05/2000 *
* *
* First, the function start_randomic() should be called to create *
* the seed (odd number): *
* seed = start_randomic(); *
* *
* Return: start_randomic -> unsigned long int *
* FRANDOM -> double in [0,1) *
********************************************************************/
unsigned rand4init(void)
{
unsigned long long y;
y = (zseed*16807LL);
zseed = (y&0x7fffffff) + (y>>31);
if (zseed&0x80000000)
zseed = (zseed&0x7fffffff) + 1;
return zseed;
}
void Init_Random(void)
{
int i;
ip=128;
ip1=ip-24;
ip2=ip-55;
ip3=ip-61;
for (i=ip3; i<ip; i++)
ira[i] = rand4init();
}
/****************************************************************
* Random Seed *
* Last Modified: 18/05/2000 *
* *
* When debugging, use always the same seed, otherwise, take it *
* as the numbers of seconds ellapsed since ?/?/1970. *
****************************************************************/
unsigned long start_randomic(void)
{
unsigned long semente;
semente = (unsigned long) time(NULL);
if (semente%2==0) ++semente; /*odd number*/ /* random seed */
#ifdef DEBUG
semente = 1206544417 /* 1189007899 */ /*2067511789*//*2067512563*//*1852686583*/;
#endif
zseed = semente;
Init_Random();
return semente;
}
/****************************************************************
* Random Seed *
* Last Modified: 01/04/2008 *
* *
* When debugging, use always the same seed, otherwise, take it *
* as the numbers of seconds ellapsed since ?/?/1970. *
****************************************************************/
void start_randomic_parallel(unsigned long semente)
{
zseed = semente;
Init_Random();
return;
}
/*****************************************************************
*** 2D Neighbours ***
*** Last Modified: 24/01/1999 ***
*** ***
*** In a 2D, with periodic boundary conditions, create the ***
*** matrixes with the neighbours of all sites. ***
*** ***
*** Input: right,left,up & down -> matrixes with neighb. sites***
*** l_size -> linear dimension (l_size x l_size) ***
*** ***
** From right to left, top to bottom, front to back. ***
*****************************************************************/
void neighbours_2d(unsigned long *right,unsigned long *left,
unsigned long *up, unsigned long *down,
unsigned long lsize)
{
unsigned long i,l2;
l2 = lsize*lsize;
for (i=0; i<l2; ++i)
{
if (i % lsize==lsize-1) *(right+i) = i-lsize+1;/* last col. */
else *(right+i) = i+1;
if (i % lsize==0) *(left+i) = i+lsize-1; /* first col.*/
else *(left+i) = i-1;
if (i<lsize) *(up+i) = l2-lsize+i; /* first row */
else *(up+i) = i - lsize;
if (i>=l2-lsize) *(down+i) = (i % lsize); /* last row */
else *(down+i) = i + lsize;
}
return;
}
/*****************************************************************
*** 2D Hexagonal Neighbours ***
*** Last Modified: 13/08/2012 ***
*** ***
*** In a 2D, with periodic boundary conditions, create the ***
*** matrixes with the neighbours of all sites. (Honeycomb) ***
*** ***
*** Input: right,left,up1,up2,down1 & down2 -> matrixes with ***
*** neighb. sites ***
*** l_size -> linear dimension (l_size x l_size) ***
*** ***
** From right to left, top to bottom, ***
*****************************************************************/
void neighbours_hexagonal_2d(unsigned long *right,unsigned long *left,
unsigned long *up1, unsigned long *down1,
unsigned long *up2, unsigned long *down2,
unsigned long lsize)
{
unsigned long i,l2;
if (lsize%2==1)
{
printf("For honeycomb lattice, LSIZE must be even.\n");
exit(8);
}
else
{
unsigned long *up,*down;
l2 = lsize*lsize;
up=(unsigned long *)calloc(l2,sizeof(unsigned long));
down=(unsigned long *)calloc(l2,sizeof(unsigned long));
if ((up!=NULL) && (down!=NULL))
{
/*Regular Square lattice*/
for (i=0; i<l2; ++i)
{
if (i % lsize==lsize-1) *(right+i) = i-lsize+1;/* last col. */
else *(right+i) = i+1;
if (i % lsize==0) *(left+i) = i+lsize-1; /* first col.*/
else *(left+i) = i-1;
if (i<lsize) *(up+i) = l2-lsize+i; /* first row */
else *(up+i) = i - lsize;
if (i>=l2-lsize) *(down+i) = (i % lsize); /* last row */
else *(down+i) = i + lsize;
}
/*Modification to hexagonal lattice*/
for (i=0; i<l2; ++i)
{
if ((i/lsize)%2==1)
{
*(up1+i)=*(up+*(left+i));
*(up2+i)=*(up+i);
*(down1+i)=*(down+*(left+i));
*(down2+i)=*(down+i);
}
else
{
*(up1+i)=*(up+i);
*(up2+i)=*(up+*(right+i));
*(down1+i)=*(down+i);
*(down2+i)=*(down+*(right+i));
}
}
free(up);
free(down);
}
else
{
fprintf(stderr,"Error in routine neighbour_hexagonal_2d\n");
exit(8);
}
}
return;
}
/*****************************************************************
*** 3D Neighbours ***
*** Last Modified: 23/03/1999 ***
*** ***
*** In a 3D, with periodic boundary conditions, create the ***
*** matrixes with the neighbours of all sites. ***
*** ***
*** Input: right,left,up,down,front, back -> neighb. sites ***
*** lsize -> linear dimension (lsize^3) ***
*****************************************************************/
void neighbours_3d(unsigned long *right,unsigned long *left,
unsigned long *up, unsigned long *down,
unsigned long *front, unsigned long *back,
unsigned long lsize)
{
unsigned long i,l2,l3;
l2 = lsize*lsize;
l3 = l2*lsize;
for (i=0; i<l3; ++i)
{
if (i % lsize==lsize-1) *(right+i) = i - lsize + 1; /* righmost plane */
else *(right+i) = i+1;
if (i % lsize==0) *(left+i) = i + lsize - 1; /* leftmost plane */
else *(left+i) = i-1;
if (i % l2 < lsize) *(up+i) = l2 - lsize + i; /* top plane */
else *(up+i) = i - lsize;
if (i % l2 >= (l2-lsize)) *(down+i) = i + lsize - l2; /* bottom plane */
else *(down+i) = i + lsize;
if (i < l2) *(front+i) = l3 - l2 + i; /* frontal plane */
else *(front+i) = i - l2;
if (i >= (l3-l2)) *(back+i) = i % l2;
else *(back+i) = i + l2; /* back plane */
}
return;
}
/*****************************************************************
*** 4D Neighbours ***
*** Last Modified: 09/05/2002 ***
*** ***
*** In a 4D, with periodic boundary conditions, create the ***
*** matrixes with the neighbours of all sites. ***
*** ***
*** Input: right,left,up,down,front,back, ***
*** front4d,back4d -> neighb. sites ***
*** lsize -> linear dimension (lsize^4) ***
*****************************************************************/
void neighbours_4d(unsigned long *right,unsigned long *left,
unsigned long *up, unsigned long *down,
unsigned long *front, unsigned long *back,
unsigned long *front4d, unsigned long *back4d,
unsigned long lsize)
{
unsigned long i,l2,l3,l4;
l2 = lsize*lsize;
l3 = l2*lsize;
l4 = l2*l2;
for (i=0; i<l4; ++i)
{
if (i % lsize==lsize-1) *(right+i) = i - lsize + 1; /* righmost plane */
else *(right+i) = i+1;
if (i % lsize==0) *(left+i) = i + lsize - 1; /* leftmost plane */
else *(left+i) = i-1;
if (i % l2 < lsize) *(up+i) = l2 - lsize + i; /* top plane */
else *(up+i) = i - lsize;
if (i % l2 >= (l2-lsize)) *(down+i) = i + lsize - l2; /* bottom plane */
else *(down+i) = i + lsize;
if (i % l3 < l2) *(front+i) = l3 - l2 + i; /* frontal plane */
else *(front+i) = i - l2;
if (i % l3 >= (l3-l2)) *(back+i) = i - (lsize-1)*l2;
else *(back+i) = i + l2; /* back plane */
if (i < l3) *(front4d+i) = l4 - l3 + i; /* 'front' cube */
else *(front4d+i) = i - l3;
if (i >=(l4-l3)) *(back4d+i) = i % l3; /* 'back' cube */
else *(back4d+i) = i + l3;
}
return;
}
/*******************************************************************************
* Instantaneous Percolation Measurements *
* Last modified: 13/06/2005 *
* *
* Remember that a cluster is identified by the last site in the chain of *
* pointers. *
*******************************************************************************/
void unionfind(unsigned long i,unsigned long j, unsigned long *lab)
{
unsigned long i1,j1;
i1 = lab[i]; /* check where i points to */
while (lab[i1] != i1) /* while it doesn't point to itself */
i1 = lab[i1];
j1 = lab[j]; /* check where j points to */
while (lab[j1] != j1) /* while it doesn't point to itself */
j1 = lab[j1];
if (lab[i1] > lab[j1]) lab[i1] = lab[j1];
else lab[j1] = lab[i1];
return;
}
/**********************************************************************
*** Hamming Distance ***
*** Last modified: 24/01/1999 ***
*** Returns the number of different sites between two vectors. ***
**********************************************************************/
int hamming_distance(int *s1, int *s2, int lsize)
{
unsigned long i,temp=0;
for (i=0; i<lsize; ++i)
if (s1[i] != s2[i]) ++temp;
return temp;
}
/*********************************************************************
*** Magnetization ***
*** Last Modified: 30/05/1999 ***
*** ***
*** Return: LSIZE*magnetization ***
*********************************************************************/
int magnetization(int *s, int lsize)
{
unsigned long i,temp=0;
for (i=0; i<lsize; ++i)
temp += s[i];
return temp;
}
int magnetization_diluted(int *s, int *n, int lsize)
{
unsigned long i, temp=0;
for (i=0; i<lsize; ++i)
temp += (s[i]*n[i]);
return temp;
}
/*********************************************************************
*** Overlap ***
*** Last Modified: 30/05/1999 ***
*** ***
*** Return: LSIZE*overlap ***
*********************************************************************/
int overlap(int *s1, int *s2, int lsize)
{
unsigned long i,temp=0;
for (i=0; i<lsize; ++i)
temp += s1[i]*s2[i];
return temp;
}
int overlap_diluted(int *s1, int *n1, int *s2, int *n2, int lsize)
{
unsigned long i,temp=0;
for (i=0; i<lsize; ++i)
temp += s1[i]*n1[i]*s2[i]*n2[i];
return temp;
}
/********************************************************************
*** Gaussian Random Number Generator ***
*** Last Modified: 18/05/2000 ***
********************************************************************/
float ngaussian(void)
{
static int iset=0;
static float gset;
float fac,r,v1,v2;
if (iset==0) {
do {
v1=2.0*FRANDOM-1.0;
v2=2.0*FRANDOM-1.0;
r=v1*v1+v2*v2;
}
while (r>=1 || r==0.0);
fac=sqrt(-2.0*log(r)/r);
gset=v1*fac;
iset=1;
return v2*fac;
}
else {
iset=0;
return gset;
}
}
/*********************************************************************
*** Create Initial Bond Configuration ***
*** Last modified: 18/05/2000 ***
*** ***
*** The (symmetric) bonds are created depending on the model ***
*** up to dimension 3: ***
*** 0. Ferromagnet: all bonds = 1 ***
*** 1. Spin Glass: same probability of being +1 and -1 ***
*** 2. Spin Glass: gaussian distributed ***
*** 3. Fully Frustrated: horizontal connections = +1 and ***
*** vertical ones alternating signals such that all minimal ***
*** loops are frustrated (LSIZE must be EVEN!) ***
*** 2D ONLY!! ***
*** ***
*** 4. Fully Frustrated (Federico) ***
*** ***
*** The first index follows the rule: ***
*** 0. Horizontal coupling; ***
*** 1. Vertical coupling; ***
*** 2. Back coupling; ***
*** ***
*** Use 3 for the FF-2D! ***
*********************************************************************/
void create_bonds(int model, int coordination, int lsize,
double *connect)
{
unsigned long i,k,ll,x,y,z; //l2;
ll = lsize;
// l2 = lsize*lsize;
switch (coordination)
{
case 2: ll *= lsize; /* 2D */
break;
case 3: ll *= lsize*lsize; /* 3D */
break;
default: fprintf(stdout,"Dimension should be either 2 or 3!\n");
exit(1);
break;
}
for (k=0; k<coordination; ++k)
for (i=0; i<ll; ++i)
switch (model){
case 0: *(connect+k*ll+i) = 1;
break;
case 1: if (FRANDOM<0.5) *(connect+k*ll+i) = 1;
else *(connect+k*ll+i) = -1;
break;
case 2: *(connect+k*ll+i) = ngaussian();
break;
case 3: *(connect+i) = 1; /* connect[0][i] */
if ((i % lsize) % 2 == 0) *(connect+ll+i) = -1;
else *(connect+ll+i) = 1;
break;
case 4: x = i % lsize; /* From Federico */
y = ((i - x) / lsize) % lsize;
z = i / (lsize*lsize);
*(connect+i) = 1;
*(connect+ll+i) = 1;
*(connect+2*ll+i) = 1;
if ((y % 2 == 0) && (z % 2 == 1))
*(connect+i) = -1;
if ((x % 2 == 0) && (z % 2 == 0))
*(connect+ll+i) = -1;
if ((x % 2 == 1) && (y % 2 == 1))
*(connect+2*ll+i) = -1;
break;
default: fprintf(stdout,"Bond set not defined!\n");
exit(1);
break;
}
return;
}
void create_bonds_a(double aa, int coordination, int lsize, double *connect)
{
unsigned long i,k,ll;
ll = lsize*lsize;
if (coordination==3) ll *= lsize;
else if (coordination!=2)
fprintf(stdout,"Dimension should be either 2 or 3!\n");
for (k=0; k<coordination; ++k)
for (i=0; i<ll; ++i)
{
if (FRANDOM<aa) *(connect+k*ll+i) = 1;
else *(connect+k*ll+i) = -1;
}
return;
}
/*********************************************************************
*** Create Initial Configuration ***
*** Last modified: 25/04/2000 ***
*** ***
*** Create at random an initial configuration for both S_i and n_i.***
*********************************************************************/
void initial_configuration(float mag, float dens, int ll,
int *s, int *n)
{
unsigned long i;
float rm;
rm = 0.5*(1.0 - mag); /* fraction of -1 sites */
if (dens!=0.0)
for (i=0; i<ll; ++i)
{
if (s != NULL)
{
if (FRANDOM<rm) s[i] = -1;
else s[i] = 1;
}
if (n != NULL)
{
if (FRANDOM<dens) n[i] = 1;
else n[i] = 0;
}
}
return;
}
/*********************************************************************
*** Measuring Performance (time) ***
*** Last modified: 24/03/1999 ***
*** ***
*** Use: ***
*** cpu_time = (clock() - begin_clock) / (double) CLOCKS_PER_SEC; ***
*** real_time = difftime(time(NULL), begin_time); ***
*********************************************************************/
void start_stopwatches(void)
{
begin_cpu_time = clock();
begin_real_time = time(NULL);
return;
}
/*********************************************************************
*** Simple Matrix Operations ***
*** Last Modified: 18/05/2000 ***
*** ***
*********************************************************************/
void copymatrix(int *s1, int *s2, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=s2[i];
return;
}
void copymatrix_f(float *s1, float *s2, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=s2[i];
return;
}
void zeromatrix_i(int *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0;
return;
}
void zeromatrix_l(long *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0;
return;
}
void zeromatrix_ul(unsigned long *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0;
return;
}
void zeromatrix_sl(signed long *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0;
return;
}
void zeromatrix_f(float *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0.0;
return;
}
void zeromatrix_d(double *s1, int lsize)
{
unsigned long i;
for (i=0; i<lsize; ++i)
s1[i]=0.0;
return;
}
/**********************************************************************
*** Create Time Table ***
*** Last Modified: 04/08/1999 ***
**********************************************************************/
void create_time_table(long *t1,long *t2,int power, int base)
{
long i,j;
t1[0] = 1;
for (i=1; i<power; i++)
t1[i] = t1[i-1]*base;
if (t2 != NULL)
{
for (i=0; i<power; i++)
for (j=0; j<power; j++)
{
*(t2+i*power+j) = t1[i] + t1[j];
if (*(t2+i*power+j) == base*t1[power-1]) *(t2+i*power+j) =base*t1[power-1];
if (*(t2+i*power+j) > base*t1[power-1]) *(t2+i*power+j) =-1;
}
}
return;
}
/*********************************************************************
*** Time Table 2 ***
*** Last Modified: 08/01/2000 ***
*** The total number of time steps and the number of measures are ***
*** specified. ***
*********************************************************************/
void create_time_table_2(long *t1, int total_time, int measures)
{
unsigned long i,k;
double temp;
t1[0] = 0;
temp = pow((double) total_time,1.0/(measures-1));
k=0;
for (i=1; i<measures; ++i)
{
t1[i] = (int) pow(temp, (double) i);
if (t1[i]<=k) t1[i]=k+1;
k = t1[i];
}
return;
}
/*********************************************************************
*** Time Table ***
*********************************************************************/
void create_time_table_fa(float *t1, float *t2, int decades)
{
unsigned long j,k=1;
signed int i=-3;
t1[0] = 0;
while (i <= decades)
{
t1[k] = pow(10,i-1);
++k;
for (j=2; j<10; ++j)
{
t1[k] = t1[k-1] + pow(10,i-1);
++k;
}
++i;
}
t1[k] = pow(10,i-1);
if (t2 != NULL)
{
for (i=0; i<decades+5; i++)
for (j=0; j<(decades+4)*9+2; j++)
{
if (i==0) *(t2+i*((decades+4)*9+2)+j) = t1[j];
else
*(t2+i*((decades+4)*9+2)+j) = t1[9*i+1] + t1[j];
/*if (*(t2+i*10+j) == base*t1[power-1]) *(t2+i*power+j) =base*t1[power-1];
if (*(t2+i*10+j) > base*t1[power-1]) *(t2+i*power+j) =-1;*/
}
}
return;
}