Modelo de ising com magnetização constante e algoritmo de Kawasaki
Ir para navegação
Ir para pesquisar
/* 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;
}