alelimafisica/Multicanonical-Methods — Issue #1: Wang-Landau and BHM
Author: alelimafisica
Created: 2012-02-26 17:09:06+00:00
State: open

#include "defs.h"
#include

#define e 2.71828182845904523536

int **NbrSpins, *HistE, **Nud, *EneSpin, *Spin;
double *LngE, *MeanMag, *Mean2Mag, *Mean4Mag, **Nmud, *MeanTau;
int L, N, E, E0, EF, Mag, oldspin, EneRange, flat, Samples, site;
double e0, ef, flatp, f, eta;
gsl_rng *r;

FILE *BHM;
FILE *WL;

void Boundaries(){
int i;

for (i=0;i= E0){
site = N * gsl_rng_uniform(r);
Actualize(site);
}
}

void Allocate(){
int k;

E0 = e0 * 2 * N;
if((E0%2)) E0++;
EF = ef * 2 * N;
if((EF%2)) EF++;
EneRange = (EF - E0)/2 + 1;

AllocMem(Spin, N, int);
AllocMem(EneSpin, N, int);
AllocMem(HistE, EneRange, int);
AllocMem(LngE, EneRange, double);
AllocMem2(NbrSpins, N, 4, int);
AllocMem2(Nud, EneRange, 5, int);
AllocMem(MeanMag, EneRange, double);
AllocMem(Mean2Mag, EneRange, double);
AllocMem(Mean4Mag, EneRange, double);
AllocMem(MeanTau, EneRange, double);
AllocMem2(Nmud, EneRange, 5, double);
}

void Actualize(int s){

int i, j, count, Eline, dE;
oldspin = Spin[s];
Spin[s] *= -1; // Atualizando o estado do spin
Mag -= 2 * oldspin; // Atualzando a magnetizacao
Eline = (E - E0)/2;
Nud[Eline][4 - EneSpin[s]]--;
E += 4 - EneSpin[s]; // Atualizando a E_total e a E_local
EneSpin[s] = 4 - EneSpin[s];
Eline = (E - E0)/2;
Nud[Eline][4 - EneSpin[s]]++;
for(i = 0; i < 4; i++){
EneSpin[NbrSpins[s][i]] += (Spin[NbrSpins[s][i]]==oldspin)?1:-1;
Nud[Eline][4 - EneSpin[NbrSpins[s][i]]]++;
}
}

int FlatTest(int n){
int i, cont;
double x, x_min, x_max;
x = 0.0;
cont = 0;

for(i = 0;i < n; i++){
if(HistE[i]){
x += (double)HistE[i];
cont++;
}
}
x /= (double)cont;
x_max = x + flatp * x;
x_min = x - flatp * x;

for(i = 0; i < n; i++){
if(HistE[i]){
if((HistE[i] > x_max) || (HistE[i] < x_min)){
flat = 0;
return(flat);
}
}
}
flat = 1;
return(flat);
}

double EvalTau(){
int i, j, count1, count2, line1, line2, collum1, collum2;
int tau;

tau = 0;
//Varredura horizontal
for(i = 0; i < L/2; i++){
count1 = count2 = 0;
line1 = i*L;
line2 = line1 + N/2;
for(j = 0; j < L; j++){
count1 += Spin[line1+j];
count2 += Spin[line2+j];
}
if((count1 * count2) > 0)
tau++;
else if((count1==0) || (count2==0));
else
tau--;
}
//Varredura vertical
for(i = 0; i < L/2; i++){
count1 = count2 = 0;
collum1 = i;
collum2 = collum1 + L/2;
for(j = 0; j < N; j += L){
count1 += Spin[collum1+j];
count2 += Spin[collum2+j];
}
if((count1 * count2) > 0)
tau++;
else if((count1==0) || (count2==0));
else
tau--;
}
return((double)tau/(double)L);
}

void EvalMeans(){

double mag;

LngE[(E-E0)/2] += log(f);
HistE[(E-E0)/2] ++;
mag = (double)(abs(Mag))/(double)N;
MeanMag[(E-E0)/2] += mag/(double)Samples;
Mean2Mag[(E-E0)/2] += (mag*mag)/(double)Samples;
Mean4Mag[(E-E0)/2] += (mag*mag*mag*mag)/(double)Samples;
MeanTau[(E-E0)/2] += EvalTau()/(double)(Samples);
}

int main(int argc, char **argv){

int i, j, dE, E1, E2;
char name[200], name1[200];
double p, fmin;

if(argc != 6){
printf("Usa: <%s>
\n Que funciona\n", argv[0]);
return(1);
}
r = gsl_rng_alloc(gsl_rng_mt19937);
L = atoi(argv[1]);
e0 = atof(argv[2]);
ef = atof(argv[3]);
Samples = atoi(argv[4]);
flatp = atof(argv[5]);
N = L * L;

Allocate();

sprintf(name, "wl_L_%d_e0_%f_ef_%f.dat", L, e0, ef);
WL = fopen(name, "w");

sprintf(name1, "Nup_Ndn_wl_L_%d_e0_%f_ef_%f.dat", L, e0, ef);
BHM = fopen(name1, "w");

Boundaries();
Initialize();
fmin = exp(0.00000001);

while (f > fmin){
for(i = 0; i < Samples; i++){
for(j = 0; j < N; j++){
site = N * gsl_rng_uniform(r);
oldspin = Spin[site];
dE = 4 - 2 * EneSpin[site];
if((E+dE <= EF) && (E+dE >= E0)){
E1 = (int)(E - E0)/2;
E2 = (int)(E + dE - E0)/2;
if ((LngE[E1] - LngE[E2]) >log( gsl_rng_uniform(r)))
Actualize(site);
}
}
EvalMeans();
}
if(FlatTest(EneRange)){
f = sqrt(f);
ZeraHist(EneRange);
}
}

for(i = 0; i < EneRange; i++){
if(LngE[i]!=0.)
fprintf(WL, "%10.12lf\t%10.12lf\t%10.12lf\t%10.12lf\t%10.12lf\t%10.12lf\t%10.12lf\t\n", (2*i + E0) / (double)(2*N) , LngE[i], MeanMag[i], Mean2Mag[i], Mean4Mag[i], MeanTau[i], (1. - 3.*(Mean4Mag[i]/(Mean2Mag[i]*Mean2Mag[i]))));
fflush(0);
}

return(0);
}